diff --git a/dumux/common/staggeredfvproblem.hh b/dumux/common/staggeredfvproblem.hh
index bed5fa94eb55d183ba88a7964b005021416e0f9a..43ddd021c35aa9b2769b9ae1119fd95518c97658 100644
--- a/dumux/common/staggeredfvproblem.hh
+++ b/dumux/common/staggeredfvproblem.hh
@@ -69,6 +69,9 @@ class StaggeredFVProblem : public FVProblem<TypeTag>
     static constexpr auto cellCenterIdx = FVGridGeometry::cellCenterIdx();
     static constexpr auto faceIdx = FVGridGeometry::faceIdx();
 
+    static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+    static constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
+
 public:
     /*!
      * \brief Constructor
@@ -174,24 +177,24 @@ public:
 
     //! Applys the initial cell center solution
     void applyInitialCellCenterSolution(SolutionVector& sol,
-                                         const SubControlVolume& scv,
-                                         const PrimaryVariables& initSol) const
+                                        const SubControlVolume& scv,
+                                        const PrimaryVariables& initSol) const
     {
-        constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
-        constexpr auto numEq = GET_PROP_TYPE(TypeTag, ModelTraits)::numEq();
-        constexpr auto offset = numEq - numEqCellCenter;
+        // while the container within the actual solution vector holds numEqCellCenter
+        // elements, we need to specify an offset to get the correct entry of the initial solution
+        static constexpr auto offset = PrimaryVariables::dimension - numEqCellCenter;
 
-        for(auto&& i : priVarIndices_(cellCenterIdx))
-            sol[cellCenterIdx][scv.dofIndex()][i] = initSol[i + offset];
+        for(int pvIdx = 0; pvIdx < numEqCellCenter; ++pvIdx)
+            sol[cellCenterIdx][scv.dofIndex()][pvIdx] = initSol[pvIdx + offset];
     }
 
     //! Applys the initial face solution
     void applyInitialFaceSolution(SolutionVector& sol,
-                                   const SubControlVolumeFace& scvf,
-                                   const PrimaryVariables& initSol) const
+                                  const SubControlVolumeFace& scvf,
+                                  const PrimaryVariables& initSol) const
     {
-        for(auto&& i : priVarIndices_(faceIdx))
-            sol[faceIdx][scvf.dofIndex()][i] = initSol[i];
+        for(int pvIdx = 0; pvIdx < numEqFace; ++pvIdx)
+            sol[faceIdx][scvf.dofIndex()][pvIdx] = initSol[pvIdx];
     }
 
 protected:
@@ -202,23 +205,6 @@ protected:
     //! \copydoc asImp_()
     const Implementation &asImp_() const
     { return *static_cast<const Implementation *>(this); }
-
-    //! Helper function that returns an iterable range of primary variable indices.
-    //! Specialization for cell center dofs.
-    static auto priVarIndices_(typename FVGridGeometry::DofTypeIndices::CellCenterIdx)
-    {
-        constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
-        return Dune::range(0, numEqCellCenter);
-    }
-
-    //! Helper function that returns an iterable range of primary variable indices.
-    //! Specialization for face dofs.
-    static auto priVarIndices_(typename FVGridGeometry::DofTypeIndices::FaceIdx)
-    {
-        constexpr auto numEqFace = GET_PROP_VALUE(TypeTag, NumEqFace);
-        return Dune::range(0, numEqFace);
-    }
-
 };
 
 } // end namespace Dumux
diff --git a/dumux/freeflow/compositional/kepsilonncmodel.hh b/dumux/freeflow/compositional/kepsilonncmodel.hh
index 2697a5e0ae067ef3b8ce1b3bac20d4ee5df3fd7a..154e295ca6383c155bddd9e7b1bdbe1c815220be 100644
--- a/dumux/freeflow/compositional/kepsilonncmodel.hh
+++ b/dumux/freeflow/compositional/kepsilonncmodel.hh
@@ -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>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, KEpsilonIndices<dimension, nComp, phaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model
diff --git a/dumux/freeflow/compositional/komegancmodel.hh b/dumux/freeflow/compositional/komegancmodel.hh
index 5259887a2f6c2640e3e990757e433a8a897721d6..efdffd868f049c9688dc8bcf68ca7a17b36cbf31 100644
--- a/dumux/freeflow/compositional/komegancmodel.hh
+++ b/dumux/freeflow/compositional/komegancmodel.hh
@@ -56,9 +56,15 @@ NEW_TYPE_TAG(KOmegaNC, INHERITS_FROM(NavierStokesNC));
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the low-Reynolds k-epsilon multi-component model
+ *
+ * \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
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
-struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
+struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -68,7 +74,7 @@ struct KOmegaNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseId
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, KOmegaIndices<dimension, nComp>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, KOmegaIndices<dimension, nComp, fluidSystemPhaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model
diff --git a/dumux/freeflow/compositional/lowrekepsilonncmodel.hh b/dumux/freeflow/compositional/lowrekepsilonncmodel.hh
index 4918178db8a539c76e3fcba1c27dd3d8f8623ad1..c11340df6fbf314ba117422f8df5beb421e97209 100644
--- a/dumux/freeflow/compositional/lowrekepsilonncmodel.hh
+++ b/dumux/freeflow/compositional/lowrekepsilonncmodel.hh
@@ -56,9 +56,15 @@ NEW_TYPE_TAG(LowReKEpsilonNC, INHERITS_FROM(NavierStokesNC));
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the low-Reynolds k-epsilon multi-component model
+ *
+ * \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
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
-struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
+struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -68,7 +74,7 @@ struct LowReKEpsilonNCModelTraits : NavierStokesNCModelTraits<dimension, nComp,
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, LowReKEpsilonIndices<dimension, nComp>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, LowReKEpsilonIndices<dimension, nComp, fluidSystemPhaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component low-Reynolds k-epsilon model
diff --git a/dumux/freeflow/compositional/navierstokesncmodel.hh b/dumux/freeflow/compositional/navierstokesncmodel.hh
index 56f25d472ea424b3e87409828a51382e1d9c4046..0c8e4056aece1e1c5b0c2e4818ecab8b4fcf848d 100644
--- a/dumux/freeflow/compositional/navierstokesncmodel.hh
+++ b/dumux/freeflow/compositional/navierstokesncmodel.hh
@@ -76,9 +76,15 @@ namespace Dumux {
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the multi-component free-flow model
+ *
+ * \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
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useM>
-struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useM>
+struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -94,7 +100,7 @@ struct NavierStokesNCModelTraits : NavierStokesModelTraits<dimension>
     static constexpr bool enableMolecularDiffusion() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, NavierStokesIndices<dimension>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, NavierStokesIndices<dimension, fluidSystemPhaseIdx>>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -127,6 +133,10 @@ private:
     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");
+
 public:
     using type = NavierStokesNCModelTraits<dim, numComponents, phaseIdx, replaceCompEqIdx, useMoles>;
 };
diff --git a/dumux/freeflow/compositional/oneeqncmodel.hh b/dumux/freeflow/compositional/oneeqncmodel.hh
index 98827f942d428ac9407ad6baa54cf48e83dfac09..a4b379e47578612b5e0aae66e9b47c3cadf9fa8e 100644
--- a/dumux/freeflow/compositional/oneeqncmodel.hh
+++ b/dumux/freeflow/compositional/oneeqncmodel.hh
@@ -56,9 +56,15 @@ NEW_TYPE_TAG(OneEqNC, INHERITS_FROM(NavierStokesNC));
 /*!
  * \ingroup FreeflowNCModel
  * \brief Traits for the one-equation multi-component model
+ *
+ * \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
  */
-template<int dimension, int nComp, int phaseIdx, int replaceCompEqIdx, bool useMoles>
-struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx, replaceCompEqIdx, useMoles>
+template<int dimension, int nComp, int fluidSystemPhaseIdx, int replaceCompEqIdx, bool useMoles>
+struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, fluidSystemPhaseIdx, replaceCompEqIdx, useMoles>
 {
     //! There are as many momentum balance equations as dimensions
     //! and as many balance equations as components.
@@ -68,7 +74,7 @@ struct OneEqNCModelTraits : NavierStokesNCModelTraits<dimension, nComp, phaseIdx
     static constexpr bool usesTurbulenceModel() { return true; }
 
     //! the indices
-    using Indices = FreeflowNCIndices<dimension, numEq(), phaseIdx, replaceCompEqIdx, OneEqIndices<dimension, nComp>>;
+    using Indices = FreeflowNCIndices<dimension, numEq(), fluidSystemPhaseIdx, replaceCompEqIdx, OneEqIndices<dimension, nComp, fluidSystemPhaseIdx>>;
 };
 
 //!< states some specifics of the isothermal multi-component one-equation model
diff --git a/dumux/freeflow/compositional/volumevariables.hh b/dumux/freeflow/compositional/volumevariables.hh
index aa76e966acb497a973f47da85d1dd80aeef1bd46..286d6f456e42784f45515fc6f798e3e1128b42bd 100644
--- a/dumux/freeflow/compositional/volumevariables.hh
+++ b/dumux/freeflow/compositional/volumevariables.hh
@@ -156,14 +156,14 @@ public:
      * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
      *        the control volume.
      */
-    Scalar pressure() const
+    Scalar pressure(int phaseIdx = fluidSystemPhaseIdx) const
     { return fluidState_.pressure(fluidSystemPhaseIdx); }
 
     /*!
      * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
      *        control volume.
      */
-    Scalar density() const
+    Scalar density(int phaseIdx = fluidSystemPhaseIdx) const
     { return fluidState_.density(fluidSystemPhaseIdx); }
 
     /*!
@@ -187,9 +187,20 @@ public:
      * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
      *        control volume.
      */
-    Scalar viscosity() const
+    Scalar viscosity(int phaseIdx = fluidSystemPhaseIdx) const
     { return fluidState_.viscosity(fluidSystemPhaseIdx); }
 
+    /*!
+     * \brief Returns the mass fraction of a component in the phase \f$\mathrm{[-]}\f$
+     *
+     * \param phaseIdx the index of the phase
+     * \param compIdx the index of the component
+     */
+    Scalar massFraction(int phaseIdx, int compIdx) const
+    {
+        return fluidState_.massFraction(fluidSystemPhaseIdx, compIdx);
+    }
+
     /*!
      * \brief Returns the mass fraction of a component in the phase \f$\mathrm{[-]}\f$
      *
@@ -200,6 +211,17 @@ public:
         return fluidState_.massFraction(fluidSystemPhaseIdx, compIdx);
     }
 
+    /*!
+     * \brief Returns the mole fraction of a component in the phase \f$\mathrm{[-]}\f$
+     *
+     * \param phaseIdx the index of the phase
+     * \param compIdx the index of the component
+     */
+    Scalar moleFraction(int phaseIdx, int compIdx) const
+    {
+        return fluidState_.moleFraction(fluidSystemPhaseIdx, compIdx);
+    }
+
     /*!
      * \brief Returns the mole fraction of a component in the phase \f$\mathrm{[-]}\f$
      *
@@ -213,7 +235,7 @@ public:
     /*!
      * \brief Returns the mass density of a given phase \f$\mathrm{[kg/m^3]}\f$
      */
-    Scalar molarDensity() const
+    Scalar molarDensity(int phaseIdx = fluidSystemPhaseIdx) const
     {
         return fluidState_.molarDensity(fluidSystemPhaseIdx);
     }
diff --git a/dumux/freeflow/navierstokes/indices.hh b/dumux/freeflow/navierstokes/indices.hh
index f77562876659d63db1cc631e4c2eddfde161c208..9bf4c76feab9c8f834cdd2cca46ec6cb5ee27995 100644
--- a/dumux/freeflow/navierstokes/indices.hh
+++ b/dumux/freeflow/navierstokes/indices.hh
@@ -31,8 +31,9 @@ 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>
+template <int dimension, int fsPhaseIdx>
 struct NavierStokesIndices
 {
     static constexpr int dimXIdx = 0; //!< Index of the x-component of a vector of size dim
@@ -63,7 +64,7 @@ struct NavierStokesIndices
     }
 
     //! The index of the fluid phase in the fluid system (for compatibility reasons)
-    static constexpr int fluidSystemPhaseIdx = 0;
+    static constexpr int fluidSystemPhaseIdx = fsPhaseIdx;
 };
 
 } // end namespace Dumux
diff --git a/dumux/freeflow/navierstokes/model.hh b/dumux/freeflow/navierstokes/model.hh
index 3d4ef755033d4ed1764c1829c7384491e7fb834a..1dd6562323f247155e760dc9d27d88bf35f97314 100644
--- a/dumux/freeflow/navierstokes/model.hh
+++ b/dumux/freeflow/navierstokes/model.hh
@@ -70,8 +70,11 @@ namespace Dumux {
 /*!
  * \ingroup NavierStokesModel
  * \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>
+template<int dimension, int fluidSystemPhaseIdx>
 struct NavierStokesModelTraits
 {
     //! The dimension of the model
@@ -100,7 +103,7 @@ struct NavierStokesModelTraits
     static constexpr bool usesTurbulenceModel() { return false; }
 
     //! the indices
-    using Indices = NavierStokesIndices<dim()>;
+    using Indices = NavierStokesIndices<dim(), fluidSystemPhaseIdx>;
 };
 
 /*!
@@ -152,9 +155,13 @@ SET_PROP(NavierStokes, ModelTraits)
 {
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
-    static constexpr int dim = GridView::dimension;
+    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>;
+    using type = NavierStokesModelTraits<dim, phaseIdx>;
 };
 
 /*!
@@ -183,9 +190,8 @@ private:
     using FST = typename GET_PROP_TYPE(TypeTag, FluidState);
     using MT = typename GET_PROP_TYPE(TypeTag, ModelTraits);
 
-    static_assert(MT::numPhases() == 1 && MT::numComponents() == 1 &&
-                  FSY::numPhases == 1 && FSY::numComponents == 1,
-                  "The Navier-Stokes model only works with a single-phase fluid system.");
+    static_assert(!FSY::isMiscible(),
+                  "The Navier-Stokes model only works with immiscible fluid systems.");
 
     using Traits = NavierStokesVolumeVariablesTraits<PV, FSY, FST, MT>;
 public:
@@ -215,8 +221,9 @@ SET_PROP(NavierStokesNI, ModelTraits)
 {
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
-    static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = NavierStokesModelTraits<dim>;
+    static constexpr auto dim = GridView::dimension;
+    static constexpr auto fluidSystemPhaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = NavierStokesModelTraits<dim, fluidSystemPhaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 381d4f202c2b0d88506d7d8494f667d4355cdea0..d3ed268b86bc865f6ec5cbc95151f59bf9524f09 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -248,7 +248,7 @@ public:
         // Handle inflow or outflow conditions.
         // Treat the staggered half-volume adjacent to the boundary as if it was on the opposite side of the boundary.
         // The respective face's outer normal vector will point in the same direction as the scvf's one.
-        if(scvf.boundary())
+        if(scvf.boundary() && problem.boundaryTypes(element, scvf).isDirichlet(Indices::pressureIdx))
             frontalFlux += inflowOutflowBoundaryFlux_(problem, element, scvf, elemVolVars, elemFaceVars);
 
         // Account for the staggered face's area. For rectangular elements, this equals the area of the scvf
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 3f2b0a71a0d697aef792f295478c4fcdee3e1be1..9b7e6fa865b3d116fe2563e2d361cd6d401b8342 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -299,7 +299,7 @@ protected:
                 residual = problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, scvf)[Indices::velocity(scvf.directionIndex())]
                                            * extrusionFactor * scvf.area();
 
-                // treat the remaining (normal) faces of the staggered control volume
+                // treat the remaining (frontal and lateral) faces of the staggered control volume
                 residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
             }
             else if(bcTypes.isSymmetry())
diff --git a/dumux/freeflow/navierstokes/volumevariables.hh b/dumux/freeflow/navierstokes/volumevariables.hh
index 9006aebf494b0a2adb30b648501f60d2d05ef943..5cd11521ce91535bf0d72cab1062524d91b4a8f8 100644
--- a/dumux/freeflow/navierstokes/volumevariables.hh
+++ b/dumux/freeflow/navierstokes/volumevariables.hh
@@ -111,14 +111,14 @@ public:
      * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
      *        the control volume.
      */
-    Scalar pressure() const
+    Scalar pressure(int phaseIdx = fluidSystemPhaseIdx) const
     { return fluidState_.pressure(fluidSystemPhaseIdx); }
 
     /*!
      * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
      *        control volume.
      */
-    Scalar density() const
+    Scalar density(int phaseIdx = fluidSystemPhaseIdx) const
     { return fluidState_.density(fluidSystemPhaseIdx); }
 
     /*!
@@ -135,7 +135,7 @@ public:
      * \brief Returns the mass density of a given phase within the
      *        control volume.
      */
-    Scalar molarDensity() const
+    Scalar molarDensity(int phaseIdx = fluidSystemPhaseIdx) const
     {
         return fluidState_.molarDensity(fluidSystemPhaseIdx);
     }
@@ -144,7 +144,7 @@ public:
      * \brief Returns the molar mass of a given phase within the
      *        control volume.
      */
-    Scalar molarMass() const
+    Scalar molarMass(int phaseIdx = fluidSystemPhaseIdx) const
     {
         return fluidState_.averageMolarMass(fluidSystemPhaseIdx);
     }
@@ -153,7 +153,7 @@ public:
      * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
      *        control volume.
      */
-    Scalar viscosity() const
+    Scalar viscosity(int phaseIdx = fluidSystemPhaseIdx) const
     { return fluidState_.viscosity(fluidSystemPhaseIdx); }
 
     /*!
diff --git a/dumux/freeflow/rans/model.hh b/dumux/freeflow/rans/model.hh
index 4723967012f991043db9b46937d0f2687f41c24f..f52023f155778da17f0e10db62efabd6cae13013 100644
--- a/dumux/freeflow/rans/model.hh
+++ b/dumux/freeflow/rans/model.hh
@@ -65,9 +65,12 @@ SET_BOOL_PROP(RANS, EnableInertiaTerms, true); //!< Explicitly force the conside
 /*!
  * \ingroup RANSModel
  * \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>
-struct RANSModelTraits : NavierStokesModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct RANSModelTraits : NavierStokesModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The model does include a turbulence model
     static constexpr bool usesTurbulenceModel() { return true; }
@@ -79,8 +82,12 @@ 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>;
+    using type = RANSModelTraits<dim, phaseIdx>;
 };
 
 //! The specific vtk output fields
@@ -105,7 +112,9 @@ SET_PROP(RANSNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = RANSModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+
+    using IsothermalTraits = RANSModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/oneeq/indices.hh b/dumux/freeflow/rans/oneeq/indices.hh
index 4ad0551ab0708f63d56db78ab1392d46861412a9..92cd6727f4c77fefd4fb5882143b4ffa6ab16b93 100644
--- a/dumux/freeflow/rans/oneeq/indices.hh
+++ b/dumux/freeflow/rans/oneeq/indices.hh
@@ -35,9 +35,10 @@ 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>
-struct OneEqIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct OneEqIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 public:
     static constexpr auto viscosityTildeEqIdx = dimension + numComponents;
diff --git a/dumux/freeflow/rans/oneeq/model.hh b/dumux/freeflow/rans/oneeq/model.hh
index 8dde555d5aa295663aacdad905a6a442be000d43..cdf9e90c7fe1c9947a25b590dc2ccd2a4258d668 100644
--- a/dumux/freeflow/rans/oneeq/model.hh
+++ b/dumux/freeflow/rans/oneeq/model.hh
@@ -94,9 +94,12 @@ namespace Properties {
 /*!
  * \ingroup OneEqModel
  * \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>
-struct OneEqModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct OneEqModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -109,7 +112,7 @@ struct OneEqModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = OneEqIndices<dim(), numComponents()>;
+    using Indices = OneEqIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -125,8 +128,9 @@ 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>;
+    using type = OneEqModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -184,7 +188,8 @@ SET_PROP(OneEqNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = OneEqModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = OneEqModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/indices.hh b/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
index 98a2fb402e6950dcd15df878a8df6d75dec41180..f690645d28c3f43bbc6ad817ff91100b2f746c92 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/indices.hh
@@ -35,9 +35,10 @@ 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>
-struct KEpsilonIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct KEpsilonIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 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 a280b4c2c0c402345c1a583af66a6c68214b54b7..ea1a23bf8b88cf759ec8fe5a7c62818fae84404a 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/model.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/model.hh
@@ -81,9 +81,12 @@ namespace Properties {
 /*!
  * \ingroup KEpsilonModel
  * \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>
-struct KEpsilonModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct KEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -96,7 +99,7 @@ struct KEpsilonModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = KEpsilonIndices<dim(), numComponents()>;
+    using Indices = KEpsilonIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -112,8 +115,9 @@ 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>;
+    using type = KEpsilonModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -171,7 +175,8 @@ SET_PROP(KEpsilonNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = KEpsilonModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = KEpsilonModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/twoeq/komega/indices.hh b/dumux/freeflow/rans/twoeq/komega/indices.hh
index bfdfaf381f98a4ecc45f9981c437c232d131719c..0bbe07a04aa9e625849d2fd3ed584dc0b76e4ad1 100644
--- a/dumux/freeflow/rans/twoeq/komega/indices.hh
+++ b/dumux/freeflow/rans/twoeq/komega/indices.hh
@@ -35,9 +35,10 @@ 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>
-struct KOmegaIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct KOmegaIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 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 eb300623f145b4a474ec94a30e12c128808b0684..41460e66d96cc037cdf9f0b16de3334210510970 100644
--- a/dumux/freeflow/rans/twoeq/komega/model.hh
+++ b/dumux/freeflow/rans/twoeq/komega/model.hh
@@ -88,9 +88,12 @@ namespace Properties {
 /*!
  *\ingroup KOmegaModel
  * \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>
-struct KOmegaModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct KOmegaModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -103,7 +106,7 @@ struct KOmegaModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! The indices
-    using Indices = KOmegaIndices<dim(), numComponents()>;
+    using Indices = KOmegaIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -119,8 +122,9 @@ 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>;
+    using type = KOmegaModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -179,7 +183,8 @@ SET_PROP(KOmegaNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = KOmegaModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = KOmegaModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh
index 0fb86e5fd230613bc0ce7afd7110549201891fd4..85c1b14091c681c91d6f187a34e66307f010cdfd 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/indices.hh
@@ -35,9 +35,10 @@ 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>
-struct LowReKEpsilonIndices : public NavierStokesIndices<dimension>
+template<int dimension, int numComponents, int fluidSystemPhaseIdx>
+struct LowReKEpsilonIndices : public NavierStokesIndices<dimension, fluidSystemPhaseIdx>
 {
 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 3bad183a27cba7a0c97c2302d4d21cf182b0f3ca..0b56c0bb70491d45d8bcc2bcb114db0aa60b2d05 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/model.hh
@@ -97,9 +97,12 @@ namespace Properties {
 /*!
  * \ingroup LowReKEpsilonModel
  * \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>
-struct LowReKEpsilonModelTraits : RANSModelTraits<dimension>
+template<int dimension, int fluidSystemPhaseIdx>
+struct LowReKEpsilonModelTraits : RANSModelTraits<dimension, fluidSystemPhaseIdx>
 {
     //! The dimension of the model
     static constexpr int dim() { return dimension; }
@@ -112,7 +115,7 @@ struct LowReKEpsilonModelTraits : RANSModelTraits<dimension>
     static constexpr int numComponents() { return 1; }
 
     //! the indices
-    using Indices = LowReKEpsilonIndices<dim(), numComponents()>;
+    using Indices = LowReKEpsilonIndices<dim(), numComponents(), fluidSystemPhaseIdx>;
 };
 
 ///////////////////////////////////////////////////////////////////////////
@@ -128,8 +131,9 @@ 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>;
+    using type = LowReKEpsilonModelTraits<dim, phaseIdx>;
 };
 
 //! The flux variables
@@ -187,7 +191,8 @@ SET_PROP(LowReKEpsilonNI, ModelTraits)
 private:
     using GridView = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::GridView;
     static constexpr int dim = GridView::dimension;
-    using IsothermalTraits = LowReKEpsilonModelTraits<dim>;
+    static constexpr int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+    using IsothermalTraits = LowReKEpsilonModelTraits<dim, phaseIdx>;
 public:
     using type = FreeflowNIModelTraits<IsothermalTraits>;
 };
diff --git a/dumux/freeflow/volumevariables.hh b/dumux/freeflow/volumevariables.hh
index 7fe634f3e705e1aed0103aa7f54d351a0c76529c..91937670f2b6a17c0bffb83cdb3e4c405a772092 100644
--- a/dumux/freeflow/volumevariables.hh
+++ b/dumux/freeflow/volumevariables.hh
@@ -171,14 +171,14 @@ public:
      * \brief Returns the total internal energy of a phase in the
      *        sub-control volume.
      */
-    Scalar internalEnergy() const
+    Scalar internalEnergy(int phaseIdx = fluidSystemPhaseIdx) const
     { return ParentType::asImp_().fluidState().internalEnergy(fluidSystemPhaseIdx); }
 
     /*!
      * \brief Returns the total enthalpy of a phase in the sub-control
      *        volume.
      */
-    Scalar enthalpy() const
+    Scalar enthalpy(int phaseIdx = fluidSystemPhaseIdx) const
     { return ParentType::asImp_().fluidState().enthalpy(fluidSystemPhaseIdx); }
 
     /*!