diff --git a/dumux/freeflow/compositional/indices.hh b/dumux/freeflow/compositional/indices.hh
deleted file mode 100644
index 0df039d4c0373d18f047153ca5188a995e68bc40..0000000000000000000000000000000000000000
--- a/dumux/freeflow/compositional/indices.hh
+++ /dev/null
@@ -1,56 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup FreeflowNCModel
- * \copydoc Dumux::FreeflowNCIndices
- */
-#ifndef DUMUX_FREEFLOW_NC_INDICES_HH
-#define DUMUX_FREEFLOW_NC_INDICES_HH
-
-#include <dumux/freeflow/navierstokes/indices.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup FreeflowNCModel
- * \brief The common indices for the isothermal multi-component free-flow model.
- */
-template <int dimension, int numEquations,
-          int phaseIdx, int theReplaceCompEqIdx,
-          class FreeflowIndices>
-struct FreeflowNCIndices : public FreeflowIndices
-{
-public:
-    //! The index of the fluid phase in the fluid system
-    static constexpr int fluidSystemPhaseIdx = phaseIdx;
-
-    //! The index of the main component
-    static constexpr int mainCompIdx = fluidSystemPhaseIdx;
-
-    //! Index of the pressure has to equal the one of the main component
-    static constexpr int pressureIdx = FreeflowIndices::conti0EqIdx + mainCompIdx;
-
-    //! The index of the component whose mass balance will be replaced by the total one
-    static constexpr int replaceCompEqIdx = theReplaceCompEqIdx;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/freeflow/compositional/kepsilonncmodel.hh b/dumux/freeflow/compositional/kepsilonncmodel.hh
index 154e295ca6383c155bddd9e7b1bdbe1c815220be..5e529df9bf43c8ac95f8bce564a7a563c1a68b0a 100644
--- a/dumux/freeflow/compositional/kepsilonncmodel.hh
+++ b/dumux/freeflow/compositional/kepsilonncmodel.hh
@@ -57,8 +57,8 @@ NEW_TYPE_TAG(KEpsilonNC, INHERITS_FROM(NavierStokesNC));
  * \ingroup FreeflowNCModel
  * \brief Traits for the low-Reynolds k-epsilon multi-component model
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
-struct KEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, bool useMoles, int replaceCompEqIdx>
+struct KEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, useMoles, replaceCompEqIdx>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -68,7 +68,7 @@ struct KEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phase
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, KEpsilonIndices<dimension, nComp, phaseIdx>>;
+    using Indices = KEpsilonIndices<dimension, nComp>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model
@@ -79,11 +79,10 @@ private:
     static constexpr int dimension = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
 public:
-    using type = KEpsilonNCModelTraits<dimension, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    using type = KEpsilonNCModelTraits<dimension, numComponents, useMoles, replaceCompEqIdx>;
 };
 
 //! Set the volume variables property
@@ -126,10 +125,9 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using SinglePhaseVtkOutputFields = KEpsilonVtkOutputFields<FVGridGeometry>;
 public:
-    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 //////////////////////////////////////////////////////////////////////////
@@ -147,10 +145,9 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-    using IsothermalModelTraits = KEpsilonNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    using IsothermalModelTraits = KEpsilonNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalModelTraits>;
 };
@@ -195,11 +192,10 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = KEpsilonVtkOutputFields<FVGridGeometry>;
     using NonIsothermalFields = FreeflowNonIsothermalVtkOutputFields<BaseVtkOutputFields, ModelTraits>;
 public:
-    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 // \}
diff --git a/dumux/freeflow/compositional/komegancmodel.hh b/dumux/freeflow/compositional/komegancmodel.hh
index efdffd868f049c9688dc8bcf68ca7a17b36cbf31..83e340962f5a9f20a148ecc408dde6d58510ed6f 100644
--- a/dumux/freeflow/compositional/komegancmodel.hh
+++ b/dumux/freeflow/compositional/komegancmodel.hh
@@ -59,12 +59,11 @@ NEW_TYPE_TAG(KOmegaNC, INHERITS_FROM(NavierStokesNC));
  *
  * \tparam dimension The dimension of the problem
  * \tparam nComp The number of components to be considered
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
- * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  * \tparam useM Use molar or mass balances
+ * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  */
-template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
-struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, bool useMoles, int replaceCompEqIdx>
+struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, useMoles, replaceCompEqIdx>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -74,7 +73,7 @@ struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSy
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, KOmegaIndices<dimension, nComp, fluidSystemPhaseIdx>>;
+    using Indices = KOmegaIndices<dimension, nComp>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model
@@ -85,11 +84,10 @@ private:
     static constexpr int dimension = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
 public:
-    using type = KOmegaNCModelTraits<dimension, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    using type = KOmegaNCModelTraits<dimension, numComponents, useMoles, replaceCompEqIdx>;
 };
 
 //! Set the volume variables property
@@ -132,10 +130,9 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using SinglePhaseVtkOutputFields = KOmegaVtkOutputFields<FVGridGeometry>;
 public:
-    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 //////////////////////////////////////////////////////////////////////////
@@ -153,10 +150,9 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-    using IsothermalModelTraits = KOmegaNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    using IsothermalModelTraits = KOmegaNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalModelTraits>;
 };
@@ -201,11 +197,10 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = KOmegaVtkOutputFields<FVGridGeometry>;
     using NonIsothermalFields = FreeflowNonIsothermalVtkOutputFields<BaseVtkOutputFields, ModelTraits>;
 public:
-    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 // \}
diff --git a/dumux/freeflow/compositional/lowrekepsilonncmodel.hh b/dumux/freeflow/compositional/lowrekepsilonncmodel.hh
index c11340df6fbf314ba117422f8df5beb421e97209..9b110ee5b25b1468ab020df3118b0ed1eef65474 100644
--- a/dumux/freeflow/compositional/lowrekepsilonncmodel.hh
+++ b/dumux/freeflow/compositional/lowrekepsilonncmodel.hh
@@ -59,12 +59,11 @@ NEW_TYPE_TAG(LowReKEpsilonNC, INHERITS_FROM(NavierStokesNC));
  *
  * \tparam dimension The dimension of the problem
  * \tparam nComp The number of components to be considered
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
- * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  * \tparam useM Use molar or mass balances
+ * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  */
-template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
-struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, bool useMoles, int replaceCompEqIdx>
+struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, useMoles, replaceCompEqIdx>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -74,7 +73,7 @@ struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp,
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, LowReKEpsilonIndices<dimension, nComp, fluidSystemPhaseIdx>>;
+    using Indices = LowReKEpsilonIndices<dimension, nComp>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model
@@ -85,11 +84,10 @@ private:
     static constexpr int dimension = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
 public:
-    using type = LowReKEpsilonNCModelTraits<dimension, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    using type = LowReKEpsilonNCModelTraits<dimension, numComponents, useMoles, replaceCompEqIdx>;
 };
 
 //! Set the volume variables property
@@ -132,10 +130,9 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using SinglePhaseVtkOutputFields = LowReKEpsilonVtkOutputFields<FVGridGeometry>;
 public:
-    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 //////////////////////////////////////////////////////////////////////////
@@ -153,10 +150,9 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-    using IsothermalModelTraits = LowReKEpsilonNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    using IsothermalModelTraits = LowReKEpsilonNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalModelTraits>;
 };
@@ -201,11 +197,10 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = LowReKEpsilonVtkOutputFields<FVGridGeometry>;
     using NonIsothermalFields = FreeflowNonIsothermalVtkOutputFields<BaseVtkOutputFields, ModelTraits>;
 public:
-    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 // \}
diff --git a/dumux/freeflow/compositional/navierstokesncmodel.hh b/dumux/freeflow/compositional/navierstokesncmodel.hh
index 0c8e4056aece1e1c5b0c2e4818ecab8b4fcf848d..4f43143fdc3b4f8b9a1553b9c8a73408304e846b 100644
--- a/dumux/freeflow/compositional/navierstokesncmodel.hh
+++ b/dumux/freeflow/compositional/navierstokesncmodel.hh
@@ -60,7 +60,6 @@
 #include <dumux/discretization/fourierslaw.hh>
 
 #include "volumevariables.hh"
-#include "indices.hh"
 #include "localresidual.hh"
 #include "fluxvariables.hh"
 #include "vtkoutputfields.hh"
@@ -79,12 +78,11 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam nComp The number of components to be considered
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
- * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  * \tparam useM Use molar or mass balances
+ * \tparam repCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  */
-template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useM>
-struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension, fluidSystemPhaseIdx>
+template<int dimension, int nComp, bool useM, int repCompEqIdx = nComp>
+struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -99,8 +97,11 @@ struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension, fluidSyste
     //! The one-phase model has no molecular diffusion
     static constexpr bool enableMolecularDiffusion() { return true; }
 
+    //! Index of of a component balance eq. to be replaced by a total mass/mole balance
+    static constexpr int replaceCompEqIdx() { return repCompEqIdx; }
+
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, NavierStokesIndices<dimension, fluidSystemPhaseIdx>>;
+    using Indices = NavierStokesIndices<dimension>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -130,18 +131,13 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-
-    static_assert(phaseIdx >= 0 && phaseIdx < FluidSystem::numPhases,
-                  "PhaseIdx must be non-negative and smaller than the number of phases");
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
 
 public:
-    using type = NavierStokesNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    using type = NavierStokesNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 };
 
-SET_INT_PROP(NavierStokesNC, PhaseIdx, 0); //!< Defines the phaseIdx
 SET_BOOL_PROP(NavierStokesNC, UseMoles, false); //!< Defines whether molar (true) or mass (false) density is used
 SET_INT_PROP(NavierStokesNC, ReplaceCompEqIdx, 0); //<! Set the ReplaceCompEqIdx to 0 by default
 SET_BOOL_PROP(NavierStokesNC, EnableInertiaTerms, true); //!< Consider inertia terms by default
@@ -177,10 +173,9 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = NavierStokesVtkOutputFields<FVGridGeometry>;
 public:
-    using type = FreeflowNCVtkOutputFields<BaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<BaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 /*!
@@ -213,10 +208,9 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-    using IsothermalModelTraits = NavierStokesNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    using IsothermalModelTraits = NavierStokesNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalModelTraits>;
 };
@@ -228,11 +222,10 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = NavierStokesVtkOutputFields<FVGridGeometry>;
     using NonIsothermalFields = FreeflowNonIsothermalVtkOutputFields<BaseVtkOutputFields, ModelTraits>;
 public:
-    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 //! Use Fourier's Law as default heat conduction type
diff --git a/dumux/freeflow/compositional/oneeqncmodel.hh b/dumux/freeflow/compositional/oneeqncmodel.hh
index a4b379e47578612b5e0aae66e9b47c3cadf9fa8e..c9092e38885ed75986fba53e239d7b2ac6734bc4 100644
--- a/dumux/freeflow/compositional/oneeqncmodel.hh
+++ b/dumux/freeflow/compositional/oneeqncmodel.hh
@@ -59,12 +59,11 @@ NEW_TYPE_TAG(OneEqNC, INHERITS_FROM(NavierStokesNC));
  *
  * \tparam dimension The dimension of the problem
  * \tparam nComp The number of components to be considered
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
- * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  * \tparam useM Use molar or mass balances
+ * \tparam replaceCompEqIdx The index of the component balance equation that should be replaced by a total mass/mole balance
  */
-template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
-struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, bool useMoles, int replaceCompEqIdx>
+struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, useMoles, replaceCompEqIdx>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -74,7 +73,7 @@ struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSys
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, OneEqIndices<dimension, nComp, fluidSystemPhaseIdx>>;
+    using Indices = OneEqIndices<dimension, nComp>;
 };
 
 //!< states some specifics of the isothermal multi-component one-equation model
@@ -85,11 +84,10 @@ private:
     static constexpr int dimension = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
 public:
-    using type = OneEqNCModelTraits<dimension, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    using type = OneEqNCModelTraits<dimension, numComponents, useMoles, replaceCompEqIdx>;
 };
 
 //! Set the volume variables property
@@ -132,10 +130,9 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using SinglePhaseVtkOutputFields = OneEqVtkOutputFields<FVGridGeometry>;
 public:
-    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<SinglePhaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 //////////////////////////////////////////////////////////////////////////
@@ -153,10 +150,9 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-    using IsothermalModelTraits = OneEqNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    using IsothermalModelTraits = OneEqNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalModelTraits>;
 };
@@ -201,11 +197,10 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = OneEqVtkOutputFields<FVGridGeometry>;
     using NonIsothermalFields = FreeflowNonIsothermalVtkOutputFields<BaseVtkOutputFields, ModelTraits>;
 public:
-    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 // \}
diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh
index 3cb661ddf574d72a7b2075cfbf3d28a7f4994cc6..4bfaea185451a13f7ad2ed4813887ecef34b39cc 100644
--- a/dumux/freeflow/compositional/staggered/fluxvariables.hh
+++ b/dumux/freeflow/compositional/staggered/fluxvariables.hh
@@ -52,11 +52,11 @@ class FreeflowNCFluxVariablesImpl<TypeTag, DiscretizationMethod::staggered>
     using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
 public:
-    static constexpr auto numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
-    static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    static constexpr auto numComponents = ModelTraits::numComponents();
+    static constexpr bool useMoles = ModelTraits::useMoles();
     using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
 
     /*!
@@ -88,9 +88,9 @@ public:
         flux += MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf);
 
         // in case one balance is substituted by the total mass balance
-        if (Indices::replaceCompEqIdx < numComponents)
+        if (ModelTraits::replaceCompEqIdx() < numComponents)
         {
-            flux[Indices::replaceCompEqIdx] = std::accumulate(flux.begin(), flux.end(), 0.0);
+            flux[ModelTraits::replaceCompEqIdx()] = std::accumulate(flux.begin(), flux.end(), 0.0);
         }
 
         return flux;
diff --git a/dumux/freeflow/compositional/staggered/localresidual.hh b/dumux/freeflow/compositional/staggered/localresidual.hh
index e7df94c045184ce677f161fb4b5fa80acf54a569..710b90c807ae8f5b7909f3b0781e8cb48325829f 100644
--- a/dumux/freeflow/compositional/staggered/localresidual.hh
+++ b/dumux/freeflow/compositional/staggered/localresidual.hh
@@ -52,11 +52,11 @@ class FreeflowNCResidualImpl<TypeTag, DiscretizationMethod::staggered>
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
-    using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+    using Indices = typename ModelTraits::Indices;
 
     using CellCenterResidual = CellCenterPrimaryVariables;
 
-    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
     static constexpr int numComponents =ModelTraits::numComponents();
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
@@ -85,13 +85,13 @@ public:
             const Scalar massOrMoleFraction = useMoles ? volVars.moleFraction(compIdx) : volVars.massFraction(compIdx);
             const Scalar s =  density * massOrMoleFraction;
 
-            if (eqIdx != Indices::replaceCompEqIdx)
+            if (eqIdx != ModelTraits::replaceCompEqIdx())
                 storage[eqIdx] += s;
         }
 
         // in case one balance is substituted by the total mass balance
-        if(Indices::replaceCompEqIdx < numComponents)
-            storage[Indices::replaceCompEqIdx] = density;
+        if(ModelTraits::replaceCompEqIdx() < numComponents)
+            storage[ModelTraits::replaceCompEqIdx()] = density;
 
         EnergyLocalResidual::fluidPhaseStorage(storage, volVars);
 
diff --git a/dumux/freeflow/compositional/vtkoutputfields.hh b/dumux/freeflow/compositional/vtkoutputfields.hh
index e0bbce332913632741dea8f53404eeee71e05a95..8c311dc65581f92fa46724fa9329ac23673aedcd 100644
--- a/dumux/freeflow/compositional/vtkoutputfields.hh
+++ b/dumux/freeflow/compositional/vtkoutputfields.hh
@@ -33,7 +33,7 @@ namespace Dumux
  * \ingroup FreeflowNCModel
  * \brief Adds vtk output fields specific to the FreeflowNC model
  */
-template<class BaseVtkOutputFields, class ModelTraits, class FVGridGeometry, class FluidSystem, int phaseIdx>
+template<class BaseVtkOutputFields, class ModelTraits, class FVGridGeometry, class FluidSystem>
 class FreeflowNCVtkOutputFields
 {
 
@@ -52,15 +52,16 @@ public:
     {
         for (int j = 0; j < FluidSystem::numComponents; ++j)
         {
-            vtk.addVolumeVariable([j](const auto& v){ return v.massFraction(j); }, "X^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx));
-            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.massFraction(j); }, "X^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(0));
+            vtk.addVolumeVariable([j](const auto& v){ return v.moleFraction(j); }, "x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(0));
+
+            if (j != FluidSystem::getMainComponent(0))
             {
-                vtk.addVolumeVariable([j](const auto& v){ return v.diffusionCoefficient(0, 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(0));
 
-                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())
+                    vtk.addVolumeVariable([j](const auto& v){ return v.effectiveDiffusivity(0, j) - v.diffusionCoefficient(0, j); }, "D_t");
             }
         }
     }
diff --git a/dumux/freeflow/compositional/zeroeqncmodel.hh b/dumux/freeflow/compositional/zeroeqncmodel.hh
index 61f4cc680e8e5aadb8432f47ddc7ee354cbbd29c..305aefbf0b62073c4281aba0a36ed718f498b122 100644
--- a/dumux/freeflow/compositional/zeroeqncmodel.hh
+++ b/dumux/freeflow/compositional/zeroeqncmodel.hh
@@ -57,8 +57,8 @@ NEW_TYPE_TAG(ZeroEqNC, INHERITS_FROM(NavierStokesNC));
  * \ingroup ZeroEqModel
  * \brief Traits for the Reynolds-averaged Navier-Stokes 0-Eq. model
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useM>
-struct ZeroEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useM>
+template<int dimension, int nComp, bool useM, int replaceCompEqIdx>
+struct ZeroEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, useM, replaceCompEqIdx>
 {
     //! The model does include a turbulence model
     static constexpr bool usesTurbulenceModel() { return true; }
@@ -72,11 +72,10 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
 public:
-    using type = ZeroEqNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    using type = ZeroEqNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 };
 
 //! Set the volume variables property
@@ -101,10 +100,9 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = RANSVtkOutputFields<FVGridGeometry>;
 public:
-     using type = FreeflowNCVtkOutputFields<BaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+     using type = FreeflowNCVtkOutputFields<BaseVtkOutputFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 //////////////////////////////////////////////////////////////////////////
@@ -122,10 +120,9 @@ private:
     static constexpr int dim = GridView::dimension;
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
     static constexpr int numComponents = FluidSystem::numComponents;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
     static constexpr bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-    using IsothermalModelTraits = ZeroEqNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
+    static constexpr int replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx);
+    using IsothermalModelTraits = ZeroEqNCModelTraits<dim, numComponents, useMoles, replaceCompEqIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalModelTraits>;
 };
@@ -152,11 +149,10 @@ private:
     using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
     using BaseVtkOutputFields = RANSVtkOutputFields<FVGridGeometry>;
     using NonIsothermalFields = FreeflowNonIsothermalVtkOutputFields<BaseVtkOutputFields, ModelTraits>;
 public:
-    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem, phaseIdx>;
+    using type = FreeflowNCVtkOutputFields<NonIsothermalFields, ModelTraits, FVGridGeometry, FluidSystem>;
 };
 
 // \}
diff --git a/dumux/freeflow/navierstokes/indices.hh b/dumux/freeflow/navierstokes/indices.hh
index 9bf4c76feab9c8f834cdd2cca46ec6cb5ee27995..13cd017e69c08493b2edccad20f9157a44b500d9 100644
--- a/dumux/freeflow/navierstokes/indices.hh
+++ b/dumux/freeflow/navierstokes/indices.hh
@@ -31,9 +31,8 @@ namespace Dumux {
  * \brief The common indices for the isothermal Navier-Stokes model.
  *
  * \tparam dimension The dimension of the problem
- * \tparam fsPhaseIdx The the index of the phase used for the fluid system
  */
-template <int dimension, int fsPhaseIdx>
+template <int dimension>
 struct NavierStokesIndices
 {
     static constexpr int dimXIdx = 0; //!< Index of the x-component of a vector of size dim
@@ -62,9 +61,6 @@ struct NavierStokesIndices
     {
         return dirIdx;
     }
-
-    //! The index of the fluid phase in the fluid system (for compatibility reasons)
-    static constexpr int fluidSystemPhaseIdx = fsPhaseIdx;
 };
 
 } // end namespace Dumux
diff --git a/dumux/freeflow/navierstokes/model.hh b/dumux/freeflow/navierstokes/model.hh
index 1dd6562323f247155e760dc9d27d88bf35f97314..7cfce0fd9f7b0dede4b86f5ba82bf19ea494e054 100644
--- a/dumux/freeflow/navierstokes/model.hh
+++ b/dumux/freeflow/navierstokes/model.hh
@@ -72,9 +72,8 @@ namespace Dumux {
  * \brief Traits for the Navier-Stokes model
  *
  * \tparam dimension The dimension of the problem
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int fluidSystemPhaseIdx>
+template<int dimension>
 struct NavierStokesModelTraits
 {
     //! The dimension of the model
@@ -103,7 +102,7 @@ struct NavierStokesModelTraits
     static constexpr bool usesTurbulenceModel() { return false; }
 
     //! the indices
-    using Indices = NavierStokesIndices<dim(), fluidSystemPhaseIdx>;
+    using Indices = NavierStokesIndices<dim()>;
 };
 
 /*!
@@ -156,12 +155,8 @@ SET_PROP(NavierStokes, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr auto dim = GridView::dimension;
-    static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-
-    static_assert(phaseIdx >= 0 && phaseIdx < GET_PROP_TYPE(TypeTag, FluidSystem)::numPhases,
-                  "PhaseIdx must be non-negative and smaller than the number of phases");
 public:
-    using type = NavierStokesModelTraits<dim, phaseIdx>;
+    using type = NavierStokesModelTraits<dim>;
 };
 
 /*!
@@ -222,8 +217,7 @@ SET_PROP(NavierStokesNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr auto dim = GridView::dimension;
-    static constexpr auto fluidSystemPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    using IsothermalTraits = NavierStokesModelTraits<dim, fluidSystemPhaseIdx>;
+    using IsothermalTraits = NavierStokesModelTraits<dim>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/model.hh b/dumux/freeflow/rans/model.hh
index f52023f155778da17f0e10db62efabd6cae13013..3532e1320a7e0e8b12404c8ec0e6193a3bc36593 100644
--- a/dumux/freeflow/rans/model.hh
+++ b/dumux/freeflow/rans/model.hh
@@ -67,10 +67,9 @@ SET_BOOL_PROP(RANS, EnableInertiaTerms, true); //!< Explicitly force the conside
  * \brief Traits for the Reynolds-averaged Navier-Stokes model
  *
  * \tparam dimension The dimension of the problem
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int fluidSystemPhaseIdx>
-struct RANSModelTraits : NavierStokesModelTraits<dimension, fluidSystemPhaseIdx>
+template<int dimension>
+struct RANSModelTraits : NavierStokesModelTraits<dimension>
 {
     //! The model does include a turbulence model
     static constexpr bool usesTurbulenceModel() { return true; }
@@ -82,12 +81,8 @@ SET_PROP(RANS, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-
-    static_assert(phaseIdx >= 0 && phaseIdx < GET_PROP_TYPE(TypeTag, FluidSystem)::numPhases,
-                  "PhaseIdx must be non-negative and smaller than the number of phases");
 public:
-    using type = RANSModelTraits<dim, phaseIdx>;
+    using type = RANSModelTraits<dim>;
 };
 
 //! The specific vtk output fields
@@ -112,9 +107,8 @@ SET_PROP(RANSNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 
-    using IsothermalTraits = RANSModelTraits<dim, phaseIdx>;
+    using IsothermalTraits = RANSModelTraits<dim>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/oneeq/indices.hh b/dumux/freeflow/rans/oneeq/indices.hh
index 92cd6727f4c77fefd4fb5882143b4ffa6ab16b93..4ad0551ab0708f63d56db78ab1392d46861412a9 100644
--- a/dumux/freeflow/rans/oneeq/indices.hh
+++ b/dumux/freeflow/rans/oneeq/indices.hh
@@ -35,10 +35,9 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents, int fluidSystemPhaseIdx>
-struct OneEqIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
+template<int dimension, int numComponents>
+struct OneEqIndices : public NavierStokesIndices<dimension>
 {
 public:
     static constexpr auto viscosityTildeEqIdx = dimension + numComponents;
diff --git a/dumux/freeflow/rans/oneeq/model.hh b/dumux/freeflow/rans/oneeq/model.hh
index cdf9e90c7fe1c9947a25b590dc2ccd2a4258d668..1b03724aae2d2405524f06ed1f9ef38208e1abcf 100644
--- a/dumux/freeflow/rans/oneeq/model.hh
+++ b/dumux/freeflow/rans/oneeq/model.hh
@@ -96,10 +96,9 @@ namespace Properties {
  * \brief Traits for the Spalart-Allmaras model
  *
  * \tparam dimension The dimension of the problem
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int fluidSystemPhaseIdx>
-struct OneEqModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
+template<int dimension>
+struct OneEqModelTraits : RANSModelTraits<dimension>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -112,7 +111,7 @@ struct OneEqModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = OneEqIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
+    using Indices = OneEqIndices<dim(), numComponents()>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -128,9 +127,8 @@ SET_PROP(OneEq, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = OneEqModelTraits<dim, phaseIdx>;
+    using type = OneEqModelTraits<dim>;
 };
 
 //! The flux variables
@@ -188,8 +186,7 @@ SET_PROP(OneEqNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    using IsothermalTraits = OneEqModelTraits<dim, phaseIdx>;
+    using IsothermalTraits = OneEqModelTraits<dim>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/indices.hh b/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
index f690645d28c3f43bbc6ad817ff91100b2f746c92..98a2fb402e6950dcd15df878a8df6d75dec41180 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
@@ -35,10 +35,9 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents, int fluidSystemPhaseIdx>
-struct KEpsilonIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
+template<int dimension, int numComponents>
+struct KEpsilonIndices : public NavierStokesIndices<dimension>
 {
 public:
     static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/model.hh b/dumux/freeflow/rans/twoeq/kepsilon/model.hh
index ea1a23bf8b88cf759ec8fe5a7c62818fae84404a..8a0a07a17af84e0b3eec05206a7de567943f4c4d 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/model.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/model.hh
@@ -83,10 +83,9 @@ namespace Properties {
  * \brief Traits for the k-epsilon model
  *
  * \tparam dimension The dimension of the problem
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int fluidSystemPhaseIdx>
-struct KEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
+template<int dimension>
+struct KEpsilonModelTraits : RANSModelTraits<dimension>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -99,7 +98,7 @@ struct KEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = KEpsilonIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
+    using Indices = KEpsilonIndices<dim(), numComponents()>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -115,9 +114,8 @@ SET_PROP(KEpsilon, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = KEpsilonModelTraits<dim, phaseIdx>;
+    using type = KEpsilonModelTraits<dim>;
 };
 
 //! The flux variables
@@ -175,8 +173,7 @@ SET_PROP(KEpsilonNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    using IsothermalTraits = KEpsilonModelTraits<dim, phaseIdx>;
+    using IsothermalTraits = KEpsilonModelTraits<dim>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 91fe3d55be0652977dee0776958d07f1450ea3ba..c09cbc9cf9b1f711e45f99f5d647be8dde5b43c8 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -387,7 +387,7 @@ public:
         // component mass fluxes
         for (int compIdx = 0; compIdx < ModelTraits::numComponents(); ++compIdx)
         {
-            if (Indices::replaceCompEqIdx == compIdx)
+            if (ModelTraits::replaceCompEqIdx() == compIdx)
                 continue;
 
             Scalar schmidtNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity()
@@ -405,9 +405,9 @@ public:
                     + pFunction(schmidtNumber, asImp_().turbulentSchmidtNumber()));
         }
 
-        if (Indices::replaceCompEqIdx < ModelTraits::numComponents())
+        if (ModelTraits::replaceCompEqIdx() < ModelTraits::numComponents())
         {
-            wallFunctionFlux[Indices::replaceCompEqIdx] =
+            wallFunctionFlux[ModelTraits::replaceCompEqIdx()] =
                 -std::accumulate(wallFunctionFlux.begin(), wallFunctionFlux.end(), 0.0);
         }
 
diff --git a/dumux/freeflow/rans/twoeq/komega/indices.hh b/dumux/freeflow/rans/twoeq/komega/indices.hh
index 0bbe07a04aa9e625849d2fd3ed584dc0b76e4ad1..bfdfaf381f98a4ecc45f9981c437c232d131719c 100644
--- a/dumux/freeflow/rans/twoeq/komega/indices.hh
+++ b/dumux/freeflow/rans/twoeq/komega/indices.hh
@@ -35,10 +35,9 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents, int fluidSystemPhaseIdx>
-struct KOmegaIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
+template<int dimension, int numComponents>
+struct KOmegaIndices : public NavierStokesIndices<dimension>
 {
 public:
     static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;
diff --git a/dumux/freeflow/rans/twoeq/komega/model.hh b/dumux/freeflow/rans/twoeq/komega/model.hh
index 41460e66d96cc037cdf9f0b16de3334210510970..14f7055817833fac6ebdcf6816d1d93ec85f4708 100644
--- a/dumux/freeflow/rans/twoeq/komega/model.hh
+++ b/dumux/freeflow/rans/twoeq/komega/model.hh
@@ -90,10 +90,9 @@ namespace Properties {
  * \brief Traits for the k-omega model
  *
  * \tparam dimension The dimension of the problem
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int fluidSystemPhaseIdx>
-struct KOmegaModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
+template<int dimension>
+struct KOmegaModelTraits : RANSModelTraits<dimension>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -106,7 +105,7 @@ struct KOmegaModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
     static constexpr int numComponents() { return 1; }
 
     //! The indices
-    using Indices = KOmegaIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
+    using Indices = KOmegaIndices<dim(), numComponents()>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -122,9 +121,8 @@ SET_PROP(KOmega, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = KOmegaModelTraits<dim, phaseIdx>;
+    using type = KOmegaModelTraits<dim>;
 };
 
 //! The flux variables
@@ -183,8 +181,7 @@ SET_PROP(KOmegaNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    using IsothermalTraits = KOmegaModelTraits<dim, phaseIdx>;
+    using IsothermalTraits = KOmegaModelTraits<dim>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh
index 85c1b14091c681c91d6f187a34e66307f010cdfd..0fb86e5fd230613bc0ce7afd7110549201891fd4 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh
@@ -35,10 +35,9 @@ namespace Dumux {
  *
  * \tparam dimension The dimension of the problem
  * \tparam numComponents The number of considered transported components
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int numComponents, int fluidSystemPhaseIdx>
-struct LowReKEpsilonIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
+template<int dimension, int numComponents>
+struct LowReKEpsilonIndices : public NavierStokesIndices<dimension>
 {
 public:
     static constexpr auto turbulentKineticEnergyEqIdx = dimension + numComponents;
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh
index 0b56c0bb70491d45d8bcc2bcb114db0aa60b2d05..37839774a1a2a6ad2e5cd51a680ef88ad75eaa15 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh
@@ -99,10 +99,9 @@ namespace Properties {
  * \brief Traits for the low-Reynolds k-epsilon model
  *
  * \tparam dimension The dimension of the problem
- * \tparam fluidSystemPhaseIdx The the index of the phase used for the fluid system
  */
-template<int dimension, int fluidSystemPhaseIdx>
-struct LowReKEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
+template<int dimension>
+struct LowReKEpsilonModelTraits : RANSModelTraits<dimension>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -115,7 +114,7 @@ struct LowReKEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = LowReKEpsilonIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
+    using Indices = LowReKEpsilonIndices<dim(), numComponents()>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -131,9 +130,8 @@ SET_PROP(LowReKEpsilon, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
 public:
-    using type = LowReKEpsilonModelTraits<dim, phaseIdx>;
+    using type = LowReKEpsilonModelTraits<dim>;
 };
 
 //! The flux variables
@@ -191,8 +189,7 @@ SET_PROP(LowReKEpsilonNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    using IsothermalTraits = LowReKEpsilonModelTraits<dim, phaseIdx>;
+    using IsothermalTraits = LowReKEpsilonModelTraits<dim>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
index ba590e040db2ef172817ccd836834c9ab4f34ded..e2d204d4b0a4879590ca1e388f214c81cc475026 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
@@ -566,7 +566,7 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
     template<std::size_t id> using FluidSystem  = typename GET_PROP_TYPE(SubDomainTypeTag<id>, FluidSystem);
 
     static constexpr auto numComponents = GET_PROP_TYPE(SubDomainTypeTag<stokesIdx>, ModelTraits)::numComponents();
-    static constexpr auto replaceCompEqIdx = GET_PROP_TYPE(SubDomainTypeTag<stokesIdx>, ModelTraits)::Indices::replaceCompEqIdx;
+    static constexpr auto replaceCompEqIdx = GET_PROP_TYPE(SubDomainTypeTag<stokesIdx>, ModelTraits)::replaceCompEqIdx();
     static constexpr bool useMoles = GET_PROP_TYPE(SubDomainTypeTag<stokesIdx>, ModelTraits)::useMoles();
 
     static_assert(GET_PROP_TYPE(SubDomainTypeTag<darcyIdx>, ModelTraits)::numComponents() == numComponents, "Submodels must use same number of components");
diff --git a/test/freeflow/navierstokesnc/test_msfreeflow.cc b/test/freeflow/navierstokesnc/test_msfreeflow.cc
index d19b6e6960d66c0b30572e9104256c068decd14a..77c7f83a9c30030c6f52bf19be93b5e041b78c7f 100644
--- a/test/freeflow/navierstokesnc/test_msfreeflow.cc
+++ b/test/freeflow/navierstokesnc/test_msfreeflow.cc
@@ -152,7 +152,7 @@ int main(int argc, char** argv) try
 
     // intialize the vtk output module
     using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
-    StaggeredVtkOutputModule<TypeTag, GET_PROP_VALUE(TypeTag, PhaseIdx)> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+    StaggeredVtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
     VtkOutputFields::init(vtkWriter); //! Add model specific output fields
     vtkWriter.write(0.0);
 
diff --git a/test/freeflow/ransnc/test_flatplate.cc b/test/freeflow/ransnc/test_flatplate.cc
index 2767ea5126481b716f900684206053fca418bcbe..0f426e96fb7d852a60bf6a261bd1fb47c687a465 100644
--- a/test/freeflow/ransnc/test_flatplate.cc
+++ b/test/freeflow/ransnc/test_flatplate.cc
@@ -144,7 +144,7 @@ int main(int argc, char** argv) try
 
     // intialize the vtk output module
     using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
-    StaggeredVtkOutputModule<TypeTag, GET_PROP_VALUE(TypeTag, PhaseIdx)> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+    StaggeredVtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
     VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
     vtkWriter.write(0.0);