diff --git a/dumux/discretization/staggered/freeflow/fickslaw.hh b/dumux/discretization/staggered/freeflow/fickslaw.hh
index 1cd6fa434fdeb1c1e2b5d703c9f957acb29557e4..102ce6d88afbe84b9fe472bd4758bdf41b5c7828 100644
--- a/dumux/discretization/staggered/freeflow/fickslaw.hh
+++ b/dumux/discretization/staggered/freeflow/fickslaw.hh
@@ -70,8 +70,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::staggered >
     static_assert(ModelTraits::numPhases() == 1, "Only one phase allowed supported!");
 
     enum {
-        conti0EqIdx = Indices::conti0EqIdx,
-        mainCompIdx = Indices::mainCompIdx,
+        conti0EqIdx = Indices::conti0EqIdx
     };
 
 public:
@@ -97,7 +96,7 @@ public:
 
         for(int compIdx = 0; compIdx < numComponents; ++compIdx)
         {
-            if(compIdx == mainCompIdx)
+            if(compIdx == FluidSystem::getMainComponent(0))
                 continue;
 
             // get equation index
@@ -120,7 +119,7 @@ public:
         }
 
         const Scalar cumulativeFlux = std::accumulate(flux.begin(), flux.end(), 0.0);
-        flux[mainCompIdx] = - cumulativeFlux;
+        flux[FluidSystem::getMainComponent(0)] = - cumulativeFlux;
 
         // Fick's law (for binary systems) states that the net flux of moles within the bulk phase has to be zero:
         // If a given amount of molecules A travel into one direction, the same amount of molecules B have to
@@ -146,7 +145,7 @@ public:
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
 
         const Scalar insideDistance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
-        const Scalar insideD = insideVolVars.effectiveDiffusivity(compIdx);
+        const Scalar insideD = insideVolVars.effectiveDiffusivity(0, compIdx);
         const Scalar ti = calculateOmega_(insideDistance, insideD, insideVolVars.extrusionFactor());
 
         if(scvf.boundary())
@@ -156,7 +155,7 @@ public:
             const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
             const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
             const Scalar outsideDistance = (outsideScv.dofPosition() - scvf.ipGlobal()).two_norm();
-            const Scalar outsideD = outsideVolVars.effectiveDiffusivity(compIdx);
+            const Scalar outsideD = outsideVolVars.effectiveDiffusivity(0, compIdx);
             const Scalar tj = calculateOmega_(outsideDistance, outsideD, outsideVolVars.extrusionFactor());
 
             tij = scvf.area()*(ti * tj)/(ti + tj);
diff --git a/dumux/freeflow/compositional/volumevariables.hh b/dumux/freeflow/compositional/volumevariables.hh
index 286d6f456e42784f45515fc6f798e3e1128b42bd..90d576626a2babacaf83bdc9cf3af359c16f1a7c 100644
--- a/dumux/freeflow/compositional/volumevariables.hh
+++ b/dumux/freeflow/compositional/volumevariables.hh
@@ -43,7 +43,6 @@ class FreeflowNCVolumeVariables : public FreeFlowVolumeVariables< Traits, Freefl
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
 
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
     static constexpr int numComponents = Traits::ModelTraits::numComponents();
 
     static constexpr bool useMoles = Traits::ModelTraits::useMoles();
@@ -77,6 +76,7 @@ public:
 
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(fluidState_);
+
         for (unsigned int compIIdx = 0; compIIdx < numComponents; ++compIIdx)
         {
             for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
@@ -87,7 +87,7 @@ public:
                     diffCoefficient_[compIIdx][compJIdx]
                         = FluidSystem::binaryDiffusionCoefficient(fluidState_,
                                                                   paramCache,
-                                                                  fluidSystemPhaseIdx,
+                                                                  0,
                                                                   compIIdx,
                                                                   compJIdx);
                 }
@@ -106,65 +106,62 @@ public:
                                    FluidState& fluidState)
     {
         fluidState.setTemperature(ParentType::temperature(elemSol, problem, element, scv));
-        fluidState.setPressure(fluidSystemPhaseIdx, elemSol[0][Indices::pressureIdx]);
+        fluidState.setPressure(0, elemSol[0][Indices::pressureIdx]);
 
         // saturation in a single phase is always 1 and thus redundant
         // to set. But since we use the fluid state shared by the
         // immiscible multi-phase models, so we have to set it here...
-        fluidState.setSaturation(fluidSystemPhaseIdx, 1.0);
+        fluidState.setSaturation(0, 1.0);
 
         Scalar sumFracMinorComp = 0.0;
 
-        for(int compIdx = 0; compIdx < numComponents; ++compIdx)
+        for(int compIdx = 1; compIdx < numComponents; ++compIdx)
         {
-            if(compIdx == Indices::mainCompIdx)
-                continue;
-
             // temporary add 1.0 to remove spurious differences in mole fractions
             // which are below the numerical accuracy
             Scalar moleOrMassFraction = elemSol[0][Indices::conti0EqIdx+compIdx] + 1.0;
             moleOrMassFraction = moleOrMassFraction - 1.0;
             if(useMoles)
-                fluidState.setMoleFraction(fluidSystemPhaseIdx, compIdx, moleOrMassFraction);
+                fluidState.setMoleFraction(0, compIdx, moleOrMassFraction);
             else
-                fluidState.setMassFraction(fluidSystemPhaseIdx, compIdx, moleOrMassFraction);
+                fluidState.setMassFraction(0, compIdx, moleOrMassFraction);
             sumFracMinorComp += moleOrMassFraction;
         }
         if(useMoles)
-            fluidState.setMoleFraction(fluidSystemPhaseIdx, Indices::mainCompIdx, 1.0 - sumFracMinorComp);
+            fluidState.setMoleFraction(0, 0, 1.0 - sumFracMinorComp);
         else
-            fluidState.setMassFraction(fluidSystemPhaseIdx, Indices::mainCompIdx, 1.0 - sumFracMinorComp);
+            fluidState.setMassFraction(0, 0, 1.0 - sumFracMinorComp);
 
         typename FluidSystem::ParameterCache paramCache;
-        paramCache.updatePhase(fluidState, fluidSystemPhaseIdx);
+        paramCache.updateAll(fluidState);
 
-        Scalar value = FluidSystem::density(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setDensity(fluidSystemPhaseIdx, value);
+        Scalar value = FluidSystem::density(fluidState, paramCache, 0);
+        fluidState.setDensity(0, value);
 
-        value = FluidSystem::molarDensity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setMolarDensity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::molarDensity(fluidState, paramCache, 0);
+        fluidState.setMolarDensity(0, value);
 
-        value = FluidSystem::viscosity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setViscosity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::viscosity(fluidState, paramCache, 0);
+        fluidState.setViscosity(0, value);
 
         // compute and set the enthalpy
         const Scalar h = ParentType::enthalpy(fluidState, paramCache);
-        fluidState.setEnthalpy(fluidSystemPhaseIdx, h);
+        fluidState.setEnthalpy(0, h);
     }
 
     /*!
      * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
      *        the control volume.
      */
-    Scalar pressure(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.pressure(fluidSystemPhaseIdx); }
+    Scalar pressure(int phaseIdx = 0) const
+    { return fluidState_.pressure(0); }
 
     /*!
      * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
      *        control volume.
      */
-    Scalar density(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.density(fluidSystemPhaseIdx); }
+    Scalar density(int phaseIdx = 0) const
+    { return fluidState_.density(0); }
 
     /*!
      * \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
@@ -187,8 +184,8 @@ public:
      * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
      *        control volume.
      */
-    Scalar viscosity(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.viscosity(fluidSystemPhaseIdx); }
+    Scalar viscosity(int phaseIdx = 0) const
+    { return fluidState_.viscosity(0); }
 
     /*!
      * \brief Returns the mass fraction of a component in the phase \f$\mathrm{[-]}\f$
@@ -198,7 +195,7 @@ public:
      */
     Scalar massFraction(int phaseIdx, int compIdx) const
     {
-        return fluidState_.massFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.massFraction(0, compIdx);
     }
 
     /*!
@@ -208,7 +205,7 @@ public:
      */
     Scalar massFraction(int compIdx) const
     {
-        return fluidState_.massFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.massFraction(0, compIdx);
     }
 
     /*!
@@ -219,7 +216,7 @@ public:
      */
     Scalar moleFraction(int phaseIdx, int compIdx) const
     {
-        return fluidState_.moleFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.moleFraction(0, compIdx);
     }
 
     /*!
@@ -229,15 +226,15 @@ public:
      */
     Scalar moleFraction(int compIdx) const
     {
-        return fluidState_.moleFraction(fluidSystemPhaseIdx, compIdx);
+        return fluidState_.moleFraction(0, compIdx);
     }
 
     /*!
      * \brief Returns the mass density of a given phase \f$\mathrm{[kg/m^3]}\f$
      */
-    Scalar molarDensity(int phaseIdx = fluidSystemPhaseIdx) const
+    Scalar molarDensity(int phaseIdx = 0) const
     {
-        return fluidState_.molarDensity(fluidSystemPhaseIdx);
+        return fluidState_.molarDensity(0);
     }
 
     /*!
@@ -250,29 +247,30 @@ public:
         return FluidSystem::molarMass(compIdx);
     }
 
-     /*!
-     * \brief Returns the diffusion coefficient \f$\mathrm{[m^2/s]}\f$
-     *
-     * \param compIIdx the index of the component which diffusive
-     * \param compJIdx the index of the component with respect to which compIIdx diffuses
-     */
-    Scalar diffusionCoefficient(int compIIdx, int compJIdx = fluidSystemPhaseIdx) const
-    {
-        if (compIIdx == compJIdx)
-            DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for fluidSystemPhaseIdx = compIdx");
-        return diffCoefficient_[compIIdx][compJIdx];
-    }
+    /*!
+    * \brief Returns the diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+    *
+    * \param compIIdx the index of the component which diffusive
+    * \param compJIdx the index of the component with respect to which compIIdx diffuses
+    */
+   Scalar diffusionCoefficient(int compIIdx, int compJIdx = 0) const
+   {
+       if (compIIdx == compJIdx)
+           DUNE_THROW(Dune::InvalidStateException, "Diffusion coefficient called for compIIdx = compJIdx");
+       return diffCoefficient_[compIIdx][compJIdx];
+   }
+
+    /*!
+    * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+    *
+    * \param compIIdx the index of the component which diffusive
+    * \param compJIdx the index of the component with respect to which compIIdx diffuses
+    */
+   Scalar effectiveDiffusivity(int compIIdx, int compJIdx = 0) const
+   {
+       return diffusionCoefficient(compIIdx, compJIdx);
+   }
 
-     /*!
-     * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
-     *
-     * \param compIIdx the index of the component which diffusive
-     * \param compJIdx the index of the component with respect to which compIIdx diffuses
-     */
-    Scalar effectiveDiffusivity(int compIIdx, int compJIdx = fluidSystemPhaseIdx) const
-    {
-        return diffusionCoefficient(compIIdx, compJIdx);
-    }
 
     /*!
      * \brief Return the fluid state of the control volume.
@@ -282,7 +280,7 @@ public:
 
 protected:
     FluidState fluidState_;
-    std::array<std::array<Scalar, numComponents>, numComponents> diffCoefficient_;
+    std::array<std::array<Scalar, numComponents>, numComponents> diffCoefficient_ = {};
 };
 
 } // end namespace Dumux
diff --git a/dumux/freeflow/compositional/vtkoutputfields.hh b/dumux/freeflow/compositional/vtkoutputfields.hh
index 78927b11bce0ab898c747b0ecb0c1225fe18f927..e0bbce332913632741dea8f53404eeee71e05a95 100644
--- a/dumux/freeflow/compositional/vtkoutputfields.hh
+++ b/dumux/freeflow/compositional/vtkoutputfields.hh
@@ -56,12 +56,13 @@ public:
             vtk.addVolumeVariable([j](const auto& v){ return v.moleFraction(j); }, "x^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx));
             if (j != phaseIdx)
             {
-                vtk.addVolumeVariable([j](const auto& v){ return v.diffusionCoefficient(j); }, "D^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx));
+                vtk.addVolumeVariable([j](const auto& v){ return v.diffusionCoefficient(0, j); }, "D^" + FluidSystem::componentName(j) + "_" + FluidSystem::phaseName(phaseIdx));
+
+                if (ModelTraits::usesTurbulenceModel())
+                // the eddy diffusivity is recalculated for an arbitrary component which is not the phase component
+                vtk.addVolumeVariable([j](const auto& v){ return v.effectiveDiffusivity(0, j) - v.diffusionCoefficient(0, j); }, "D_t");
             }
         }
-        if (ModelTraits::usesTurbulenceModel())
-            // the eddy diffusivity is recalculated for an arbitrary component which is not the phase component
-            vtk.addVolumeVariable([](const auto& v){ return v.effectiveDiffusivity(1-phaseIdx) - v.diffusionCoefficient(1-phaseIdx); }, "D_t");
     }
 };
 
diff --git a/dumux/freeflow/navierstokes/volumevariables.hh b/dumux/freeflow/navierstokes/volumevariables.hh
index 5cd11521ce91535bf0d72cab1062524d91b4a8f8..f4391b1733fc2758fb22ef72928ce34bd54784c5 100644
--- a/dumux/freeflow/navierstokes/volumevariables.hh
+++ b/dumux/freeflow/navierstokes/volumevariables.hh
@@ -41,8 +41,6 @@ class NavierStokesVolumeVariables : public FreeFlowVolumeVariables< Traits, Navi
 
     using Scalar = typename Traits::PrimaryVariables::value_type;
 
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
-
 public:
     //! export the underlying fluid system
     using FluidSystem = typename Traits::FluidSystem;
@@ -83,43 +81,43 @@ public:
         const Scalar t = ParentType::temperature(elemSol, problem, element, scv);
         fluidState.setTemperature(t);
 
-        fluidState.setPressure(fluidSystemPhaseIdx, elemSol[0][Indices::pressureIdx]);
+        fluidState.setPressure(0, elemSol[0][Indices::pressureIdx]);
 
         // saturation in a single phase is always 1 and thus redundant
         // to set. But since we use the fluid state shared by the
         // immiscible multi-phase models, so we have to set it here...
-        fluidState.setSaturation(fluidSystemPhaseIdx, 1.0);
+        fluidState.setSaturation(0, 1.0);
 
         typename FluidSystem::ParameterCache paramCache;
-        paramCache.updatePhase(fluidState, fluidSystemPhaseIdx);
+        paramCache.updateAll(fluidState);
 
-        Scalar value = FluidSystem::density(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setDensity(fluidSystemPhaseIdx, value);
+        Scalar value = FluidSystem::density(fluidState, paramCache, 0);
+        fluidState.setDensity(0, value);
 
-        value = FluidSystem::molarDensity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setMolarDensity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::molarDensity(fluidState, paramCache, 0);
+        fluidState.setMolarDensity(0, value);
 
-        value = FluidSystem::viscosity(fluidState, paramCache, fluidSystemPhaseIdx);
-        fluidState.setViscosity(fluidSystemPhaseIdx, value);
+        value = FluidSystem::viscosity(fluidState, paramCache, 0);
+        fluidState.setViscosity(0, value);
 
         // compute and set the enthalpy
         value = ParentType::enthalpy(fluidState, paramCache);
-        fluidState.setEnthalpy(fluidSystemPhaseIdx, value);
+        fluidState.setEnthalpy(0, value);
     }
 
     /*!
      * \brief Return the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
      *        the control volume.
      */
-    Scalar pressure(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.pressure(fluidSystemPhaseIdx); }
+    Scalar pressure(int phaseIdx = 0) const
+    { return fluidState_.pressure(0); }
 
     /*!
      * \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
      *        control volume.
      */
-    Scalar density(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.density(fluidSystemPhaseIdx); }
+    Scalar density(int phaseIdx = 0) const
+    { return fluidState_.density(0); }
 
     /*!
      * \brief Return temperature \f$\mathrm{[K]}\f$ inside the sub-control volume.
@@ -135,26 +133,26 @@ public:
      * \brief Returns the mass density of a given phase within the
      *        control volume.
      */
-    Scalar molarDensity(int phaseIdx = fluidSystemPhaseIdx) const
+    Scalar molarDensity(int phaseIdx = 0) const
     {
-        return fluidState_.molarDensity(fluidSystemPhaseIdx);
+        return fluidState_.molarDensity(0);
     }
 
     /*!
      * \brief Returns the molar mass of a given phase within the
      *        control volume.
      */
-    Scalar molarMass(int phaseIdx = fluidSystemPhaseIdx) const
+    Scalar molarMass(int phaseIdx = 0) const
     {
-        return fluidState_.averageMolarMass(fluidSystemPhaseIdx);
+        return fluidState_.averageMolarMass(0);
     }
 
     /*!
      * \brief Return the dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
      *        control volume.
      */
-    Scalar viscosity(int phaseIdx = fluidSystemPhaseIdx) const
-    { return fluidState_.viscosity(fluidSystemPhaseIdx); }
+    Scalar viscosity(int phaseIdx = 0) const
+    { return fluidState_.viscosity(0); }
 
     /*!
      * \brief Return the effective dynamic viscosity \f$\mathrm{[Pa s]}\f$ of the fluid within the
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 98bf14c686aff002488b1bb1d6e8aded562d0406..91fe3d55be0652977dee0776958d07f1450ea3ba 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -391,7 +391,7 @@ public:
                 continue;
 
             Scalar schmidtNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity()
-                                   / elemVolVars[scvf.insideScvIdx()].diffusionCoefficient(compIdx);
+                                   / elemVolVars[scvf.insideScvIdx()].diffusionCoefficient(0, compIdx);
             Scalar moleToMassConversionFactor = useMoles
                                                 ? 1.0 : FluidSystem::molarMass(compIdx);
             wallFunctionFlux[compIdx] +=
diff --git a/dumux/freeflow/rans/volumevariables.hh b/dumux/freeflow/rans/volumevariables.hh
index 2920a851e8f39af2372dc396f91c245c0acd47a1..1461d93d73a4ece4399b6141b927a3edd61f8d7d 100644
--- a/dumux/freeflow/rans/volumevariables.hh
+++ b/dumux/freeflow/rans/volumevariables.hh
@@ -49,7 +49,6 @@ class RANSVolumeVariables
     using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 
     static constexpr bool enableEnergyBalance = Traits::ModelTraits::enableEnergyBalance();
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
 
 public:
 
@@ -232,13 +231,13 @@ public:
     Scalar eddyThermalConductivity() const
     { return eddyThermalConductivity_; }
 
-     /*!
-     * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
-     *
-     * \param compIIdx the index of the component which diffusive
-     * \param compJIdx the index of the component with respect to which compIIdx diffuses
-     */
-    Scalar effectiveDiffusivity(int compIIdx, int compJIdx = fluidSystemPhaseIdx) const
+    /*!
+    * \brief Returns the effective diffusion coefficient \f$\mathrm{[m^2/s]}\f$
+    *
+    * \param compIIdx the index of the component which diffusive
+    * \param compJIdx the index of the component with respect to which compIIdx diffuses
+    */
+    Scalar effectiveDiffusivity(int compIIdx, int compJIdx) const
     {
         return NavierStokesParentType::diffusionCoefficient(compIIdx, compJIdx) + eddyDiffusivity();
     }
diff --git a/dumux/freeflow/volumevariables.hh b/dumux/freeflow/volumevariables.hh
index 91937670f2b6a17c0bffb83cdb3e4c405a772092..e60b4a697c7a115a7da152dec2b94c60113e4432 100644
--- a/dumux/freeflow/volumevariables.hh
+++ b/dumux/freeflow/volumevariables.hh
@@ -137,7 +137,6 @@ class FreeFlowVolumeVariablesImplementation<Traits, Impl, true>
 {
     using ParentType = FreeFlowVolumeVariablesImplementation<Traits, Impl, false>;
     using Scalar = typename Traits::PrimaryVariables::value_type;
-    static constexpr int fluidSystemPhaseIdx = Traits::ModelTraits::Indices::fluidSystemPhaseIdx;
 
 public:
     //! export the type used for the primary variables
@@ -171,28 +170,28 @@ public:
      * \brief Returns the total internal energy of a phase in the
      *        sub-control volume.
      */
-    Scalar internalEnergy(int phaseIdx = fluidSystemPhaseIdx) const
-    { return ParentType::asImp_().fluidState().internalEnergy(fluidSystemPhaseIdx); }
+    Scalar internalEnergy(int phaseIdx = 0) const
+    { return ParentType::asImp_().fluidState().internalEnergy(0); }
 
     /*!
      * \brief Returns the total enthalpy of a phase in the sub-control
      *        volume.
      */
-    Scalar enthalpy(int phaseIdx = fluidSystemPhaseIdx) const
-    { return ParentType::asImp_().fluidState().enthalpy(fluidSystemPhaseIdx); }
+    Scalar enthalpy(int phaseIdx = 0) const
+    { return ParentType::asImp_().fluidState().enthalpy(0); }
 
     /*!
      * \brief Returns the component enthalpy \f$\mathrm{[J/(kg*K)]}\f$ in the sub-control volume.
      */
     Scalar componentEnthalpy(unsigned int compIdx) const
-    { return FluidSystem::componentEnthalpy(ParentType::asImp_().fluidState(), fluidSystemPhaseIdx, compIdx); }
+    { return FluidSystem::componentEnthalpy(ParentType::asImp_().fluidState(), 0, compIdx); }
 
     /*!
      * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
      *        of the fluid phase in the sub-control volume.
      */
     Scalar thermalConductivity() const
-    { return FluidSystem::thermalConductivity(ParentType::asImp_().fluidState(), fluidSystemPhaseIdx); }
+    { return FluidSystem::thermalConductivity(ParentType::asImp_().fluidState(), 0); }
 
     /*!
      * \brief Returns the effective thermal conductivity \f$\mathrm{[W/(m*K)]}\f$
@@ -208,7 +207,7 @@ public:
      *        in the sub-control volume.
      */
     Scalar heatCapacity() const
-    { return FluidSystem::heatCapacity(ParentType::asImp_().fluidState(), fluidSystemPhaseIdx); }
+    { return FluidSystem::heatCapacity(ParentType::asImp_().fluidState(), 0); }
 
     //! The temperature is a primary variable for non-isothermal models
     template<class ElemSol, class Problem, class Element, class Scv>
@@ -226,7 +225,7 @@ public:
     static Scalar enthalpy(const FluidState& fluidState,
                            const ParameterCache& paramCache)
     {
-        return FluidSystem::enthalpy(fluidState, paramCache, fluidSystemPhaseIdx);
+        return FluidSystem::enthalpy(fluidState, paramCache, 0);
     }
 };
 }
diff --git a/dumux/material/fluidstates/immiscible.hh b/dumux/material/fluidstates/immiscible.hh
index b1168b4026156701952f97ab53c974b08e7a756d..68b6e29cc2395dc79a167ce8912d7fe47d7e78a0 100644
--- a/dumux/material/fluidstates/immiscible.hh
+++ b/dumux/material/fluidstates/immiscible.hh
@@ -43,9 +43,6 @@ class ImmiscibleFluidState
 public:
     static constexpr int numPhases = FluidSystem::numPhases;
     static constexpr int numComponents = FluidSystem::numComponents;
-    static_assert(numPhases == numComponents,
-                  "The number of phases must be equal to the number of "
-                  "components if immiscibility is assumed!");
 
     //! export the scalar type
     using Scalar = ScalarType;
diff --git a/dumux/material/fluidsystems/1padapter.hh b/dumux/material/fluidsystems/1padapter.hh
index a7e4f40d27e04cab47948ab03107dc57f3f1d56c..3a45b138d8a47aa1f68b3313dfbdbfa138d2e12e 100644
--- a/dumux/material/fluidsystems/1padapter.hh
+++ b/dumux/material/fluidsystems/1padapter.hh
@@ -47,18 +47,15 @@ class OnePAdapter
     using ThisType = OnePAdapter<MPFluidSystem, phase>;
     using Base = BaseFluidSystem<typename MPFluidSystem::Scalar, ThisType>;
 
+    static_assert(phase < MPFluidSystem::numPhases, "Phase does not exist in multi-phase fluidsystem!");
+
     struct AdapterPolicy
     {
         using FluidSystem = MPFluidSystem;
 
         // the phase index is always zero, other phases than the chosen phase should never be called
         static int phaseIdx(int mpFluidPhaseIdx)
-        {
-            if (mpFluidPhaseIdx == phase)
-                return 0;
-            else
-                DUNE_THROW(Dune::InvalidStateException, "Only phase " << phase << " is available!");
-        }
+        { return 0; }
 
         // the main component is currently excepted to have the same index as it's phase
         // (see Fluidsystems::Base::getMainComponent for more information)
@@ -85,6 +82,8 @@ public:
 
     //! export the wrapped MultiPhaseFluidSystem type
     using MultiPhaseFluidSystem = MPFluidSystem;
+    //! the index of the phase we choose from the multi-phase fluid system
+    static constexpr int multiphaseFluidsystemPhaseIdx = phase;
 
     //! number of phases in the fluid system
     static constexpr int numPhases = 1;
@@ -100,8 +99,9 @@ public:
     /*!
      * \brief Initialize the fluid system's static parameters generically
      */
-    static void init()
-    { MultiPhaseFluidSystem::init(); }
+    template<class ...Args>
+    static void init(Args&&... args)
+    { MultiPhaseFluidSystem::init(std::forward<Args>(args)...); }
 
     /****************************************
      * Fluid phase related static parameters
@@ -220,6 +220,23 @@ public:
         return MultiPhaseFluidSystem::enthalpy(adaptFluidState(fluidState), phase);
     }
 
+    /*!
+     * \brief Returns the specific enthalpy \f$\mathrm{[J/kg]}\f$ of a component in a specific phase
+     * \param fluidState An arbitrary fluid state
+     * \param phaseIdx The index of the fluid phase to consider
+     * \param compIdx The index of the component to consider
+     *
+     */
+    template <class FluidState>
+    static Scalar componentEnthalpy(const FluidState &fluidState,
+                                    int phaseIdx,
+                                    int compIdx)
+    {
+        assert(phaseIdx == 0);
+        return MultiPhaseFluidSystem::componentEnthalpy(adaptFluidState(fluidState), phase,
+                                                        AdapterPolicy::compIdx(compIdx));
+    }
+
     using Base::viscosity;
     /*!
      * \brief The dynamic liquid viscosity \f$\mathrm{[N/m^3*s]}\f$ of the pure component.
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
index 2b314ba87cf1992b77d176ebcc395c65e2a1a924..ba590e040db2ef172817ccd836834c9ab4f34ded 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
@@ -71,6 +71,49 @@ struct StokesDarcyCouplingOptions
 
 };
 
+/*!
+ * \ingroup MultiDomain
+ * \ingroup BoundaryCoupling
+ * \ingroup StokesDarcyCoupling
+ * \brief This structs helps to check if the two sub models use the same fluidsystem.
+ * \tparam FSFF The free-flow fluidsystem
+ * \tparam FSFF The porous-medium flow fluidsystem
+ * \tparam numPhasesPM The number of phases used in the porous-medium flow model.
+ */
+template<class FSFF, class FSPM, int numPhasesPM>
+struct IsSameFluidSystem;
+
+/*!
+ * \ingroup MultiDomain
+ * \ingroup BoundaryCoupling
+ * \ingroup StokesDarcyCoupling
+ * \brief This structs helps to check if the two sub models use the same fluidsystem.
+ *        Specialization for single-phase porous-medium models.
+ * \tparam FSFF The free-flow fluidsystem
+ * \tparam FSFF The porous-medium flow fluidsystem
+ */
+template<class FSFF, class FSPM>
+struct IsSameFluidSystem<FSFF, FSPM, 1>
+{
+    static constexpr bool value = std::is_same<FSFF, FSPM>::value;
+};
+
+/*!
+ * \ingroup MultiDomain
+ * \ingroup BoundaryCoupling
+ * \ingroup StokesDarcyCoupling
+ * \brief This structs helps to check if the two sub models use the same fluidsystem.
+ *        Specialization for two-phase porous-medium models.
+ * \tparam FSFF The free-flow fluidsystem
+ * \tparam FSFF The porous-medium flow fluidsystem
+ */
+template<class FSFF, class FSPM>
+struct IsSameFluidSystem<FSFF, FSPM, 2>
+{
+    static constexpr bool value = std::is_same<typename FSFF::MultiPhaseFluidSystem, FSPM>::value;
+};
+
+
 template<class MDTraits, class CouplingManager, bool enableEnergyBalance, bool isCompositional>
 class StokesDarcyCouplingDataImplementation;
 
@@ -106,15 +149,25 @@ class StokesDarcyCouplingDataImplementationBase
     template<std::size_t id> using VolumeVariables  = typename GET_PROP_TYPE(SubDomainTypeTag<id>, GridVolumeVariables)::VolumeVariables;
     template<std::size_t id> using Problem  = typename GET_PROP_TYPE(SubDomainTypeTag<id>, Problem);
     template<std::size_t id> using FluidSystem  = typename GET_PROP_TYPE(SubDomainTypeTag<id>, FluidSystem);
+    template<std::size_t id> using ModelTraits  = typename GET_PROP_TYPE(SubDomainTypeTag<id>, ModelTraits);
 
     static constexpr auto stokesIdx = CouplingManager::stokesIdx;
     static constexpr auto darcyIdx = CouplingManager::darcyIdx;
 
+
+
+
+
+
+
     static constexpr int enableEnergyBalance = GET_PROP_TYPE(SubDomainTypeTag<stokesIdx>, ModelTraits)::enableEnergyBalance();
     static_assert(GET_PROP_TYPE(SubDomainTypeTag<darcyIdx>, ModelTraits)::enableEnergyBalance() == enableEnergyBalance,
                   "All submodels must both be either isothermal or non-isothermal");
 
-    static_assert(std::is_same<FluidSystem<stokesIdx>, FluidSystem<darcyIdx>>::value,
+    static_assert(IsSameFluidSystem<FluidSystem<stokesIdx>,
+                                    FluidSystem<darcyIdx>,
+                                    GET_PROP_TYPE(SubDomainTypeTag<darcyIdx>, ModelTraits)::numPhases()
+                                    >::value,
                   "All submodels must use the same fluid system");
 
     using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
@@ -122,10 +175,32 @@ class StokesDarcyCouplingDataImplementationBase
 public:
     StokesDarcyCouplingDataImplementationBase(const CouplingManager& couplingmanager): couplingManager_(couplingmanager) {}
 
-    /*!
-     * \brief Export the phase index used by the Stokes model.
-     */
-    static constexpr auto stokesPhaseIdx = Indices<stokesIdx>::fluidSystemPhaseIdx;
+    template<std::size_t i, bool hasAdapter = (ModelTraits<darcyIdx>::numPhases() > 1), std::enable_if_t<!hasAdapter, int> = 0>
+    static constexpr auto couplingPhaseIdx(Dune::index_constant<i>, int coupledPhaseIdx = 0)
+    { return 0; }
+
+    template<bool hasAdapter = (ModelTraits<darcyIdx>::numPhases() > 1), std::enable_if_t<hasAdapter, int> = 0>
+    static constexpr auto couplingPhaseIdx(Dune::index_constant<stokesIdx>, int coupledPhaseIdx = 0)
+    { return 0; }
+
+    template<bool hasAdapter = (ModelTraits<darcyIdx>::numPhases() > 1), std::enable_if_t<hasAdapter, int> = 0>
+    static constexpr auto couplingPhaseIdx(Dune::index_constant<darcyIdx>, int coupledPhaseIdx = 0)
+    { return FluidSystem<stokesIdx>::multiphaseFluidsystemPhaseIdx; }
+
+
+
+
+    template<std::size_t i, bool hasAdapter = (ModelTraits<darcyIdx>::numPhases() > 1), std::enable_if_t<!hasAdapter, int> = 0>
+    static constexpr auto couplingCompIdx(Dune::index_constant<i>, int coupledCompdIdx)
+    { return coupledCompdIdx; }
+
+    template<bool hasAdapter = (ModelTraits<darcyIdx>::numPhases() > 1), std::enable_if_t<hasAdapter, int> = 0>
+    static constexpr auto couplingCompIdx(Dune::index_constant<stokesIdx>, int coupledCompdIdx)
+    { return coupledCompdIdx; }
+
+    template<bool hasAdapter = (ModelTraits<darcyIdx>::numPhases() > 1), std::enable_if_t<hasAdapter, int> = 0>
+    static constexpr auto couplingCompIdx(Dune::index_constant<darcyIdx>, int coupledCompdIdx)
+    { return FluidSystem<stokesIdx>::compIdx(coupledCompdIdx); }
 
     /*!
      * \brief Returns a reference to the coupling manager.
@@ -159,9 +234,10 @@ public:
 
         Scalar momentumFlux(0.0);
         const auto& stokesContext = couplingManager_.stokesCouplingContext(scvf);
+        const auto darcyPhaseIdx = couplingPhaseIdx(darcyIdx);
 
         // - p_pm * n_pm = p_pm * n_ff
-        const Scalar darcyPressure = stokesContext.volVars.pressure(stokesPhaseIdx);
+        const Scalar darcyPressure = stokesContext.volVars.pressure(darcyPhaseIdx);
 
         if(numPhasesDarcy > 1)
         {
@@ -170,7 +246,7 @@ public:
         else // use pressure reconstruction for single phase models
         {
             const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf();
-            const Scalar pressureInterFace = scvf.directionSign() * velocity * stokesContext.volVars.viscosity(stokesPhaseIdx)/darcyPermeability(scvf) * (stokesContext.element.geometry().center() - scvf.center()).two_norm() + darcyPressure;
+            const Scalar pressureInterFace = scvf.directionSign() * velocity * stokesContext.volVars.viscosity(darcyPhaseIdx)/darcyPermeability(scvf) * (stokesContext.element.geometry().center() - scvf.center()).two_norm() + darcyPressure;
             momentumFlux = pressureInterFace;
         }
 
@@ -336,7 +412,7 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
 
 public:
     using ParentType::ParentType;
-    using ParentType::stokesPhaseIdx;
+    using ParentType::couplingPhaseIdx;
 
     /*!
      * \brief Returns the mass flux across the coupling boundary as seen from the Darcy domain.
@@ -347,7 +423,7 @@ public:
     {
         const auto& darcyContext = this->couplingManager().darcyCouplingContext(scvf);
         const Scalar velocity = darcyContext.velocity * scvf.unitOuterNormal();
-        const Scalar darcyDensity = darcyElemVolVars[scvf.insideScvIdx()].density(stokesPhaseIdx);
+        const Scalar darcyDensity = darcyElemVolVars[scvf.insideScvIdx()].density(couplingPhaseIdx(darcyIdx));
         const Scalar stokesDensity = darcyContext.volVars.density();
         const bool insideIsUpstream = velocity > 0.0;
 
@@ -365,7 +441,7 @@ public:
         const auto& stokesContext = this->couplingManager().stokesCouplingContext(scvf);
         const Scalar velocity = stokesElemFaceVars[scvf].velocitySelf();
         const Scalar stokesDensity = stokesElemVolVars[scvf.insideScvIdx()].density();
-        const Scalar darcyDensity = stokesContext.volVars.density(stokesPhaseIdx);
+        const Scalar darcyDensity = stokesContext.volVars.density(couplingPhaseIdx(darcyIdx));
         const bool insideIsUpstream = sign(velocity) == scvf.directionSign();
 
         return massFlux_(velocity * scvf.directionSign(), stokesDensity, darcyDensity, insideIsUpstream);
@@ -446,8 +522,8 @@ private:
         const auto& outsideScv = (*scvs(outsideFvGeometry).begin());
 
         // convective fluxes
-        const Scalar insideTerm = insideVolVars.density(stokesPhaseIdx) * insideVolVars.enthalpy(stokesPhaseIdx);
-        const Scalar outsideTerm = outsideVolVars.density(stokesPhaseIdx) * outsideVolVars.enthalpy(stokesPhaseIdx);
+        const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
+        const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
 
         flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
 
@@ -502,7 +578,8 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
 
 public:
     using ParentType::ParentType;
-    using ParentType::stokesPhaseIdx;
+    using ParentType::couplingPhaseIdx;
+    using ParentType::couplingCompIdx;
 
     /*!
      * \brief Returns the mass flux across the coupling boundary as seen from the Darcy domain.
@@ -619,13 +696,18 @@ protected:
 
         // treat the advective fluxes
         auto insideTerm = [&](int compIdx)
-        { return moleOrMassFraction(insideVolVars, stokesPhaseIdx, compIdx) * moleOrMassDensity(insideVolVars, stokesPhaseIdx); };
+        { return moleOrMassFraction(insideVolVars, couplingPhaseIdx(domainI), compIdx) * moleOrMassDensity(insideVolVars, couplingPhaseIdx(domainI)); };
 
         auto outsideTerm = [&](int compIdx)
-        { return moleOrMassFraction(outsideVolVars, stokesPhaseIdx, compIdx) * moleOrMassDensity(outsideVolVars, stokesPhaseIdx); };
+        { return moleOrMassFraction(outsideVolVars, couplingPhaseIdx(domainJ), compIdx) * moleOrMassDensity(outsideVolVars, couplingPhaseIdx(domainJ)); };
+
 
         for(int compIdx = 0; compIdx < numComponents; ++compIdx)
-            flux[compIdx] += this->advectiveFlux(insideTerm(compIdx), outsideTerm(compIdx), velocity, insideIsUpstream);
+        {
+            const int domainICompIdx = couplingCompIdx(domainI, compIdx);
+            const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
+            flux[domainICompIdx] += this->advectiveFlux(insideTerm(domainICompIdx), outsideTerm(domainJCompIdx), velocity, insideIsUpstream);
+        }
 
         // treat the diffusive fluxes
         const auto& insideScv = insideFvGeometry.scv(scvf.insideScvIdx());
@@ -643,7 +725,7 @@ protected:
      */
     Scalar diffusionCoefficient_(const VolumeVariables<stokesIdx>& volVars, int phaseIdx, int compIdx) const
     {
-        return volVars.effectiveDiffusivity(compIdx);
+        return volVars.effectiveDiffusivity(phaseIdx, compIdx);
     }
 
     /*!
@@ -657,6 +739,16 @@ protected:
                                                   volVars.diffusionCoefficient(phaseIdx, compIdx));
     }
 
+    Scalar getComponentEnthalpy(const VolumeVariables<stokesIdx>& volVars, int phaseIdx, int compIdx) const
+    {
+        return FluidSystem<stokesIdx>::componentEnthalpy(volVars.fluidState(), 0, compIdx);
+    }
+
+    Scalar getComponentEnthalpy(const VolumeVariables<darcyIdx>& volVars, int phaseIdx, int compIdx) const
+    {
+        return FluidSystem<darcyIdx>::componentEnthalpy(volVars.fluidState(), phaseIdx, compIdx);
+    }
+
     /*!
      * \brief Evaluate the diffusive mole/mass flux across the interface.
      */
@@ -671,37 +763,41 @@ protected:
                                         const DiffusionCoefficientAveragingType diffCoeffAvgType) const
     {
         NumEqVector diffusiveFlux(0.0);
-        const Scalar avgMolarDensity = 0.5 * volVarsI.molarDensity(stokesPhaseIdx) + 0.5 *  volVarsJ.molarDensity(stokesPhaseIdx);
+        const Scalar avgMolarDensity = 0.5 * volVarsI.molarDensity(couplingPhaseIdx(domainI)) + 0.5 *  volVarsJ.molarDensity(couplingPhaseIdx(domainJ));
         const Scalar insideDistance = this->getDistance_(scvI, scvfI);
         const Scalar outsideDistance = this->getDistance_(scvJ, scvfI);
 
-        for(int compIdx = 0; compIdx < numComponents; ++compIdx)
+        for(int compIdx = 1; compIdx < numComponents; ++compIdx)
         {
-            if(compIdx != stokesPhaseIdx)
-            {
-                const Scalar deltaMoleFrac = volVarsJ.moleFraction(stokesPhaseIdx, compIdx) - volVarsI.moleFraction(stokesPhaseIdx, compIdx);
-                const Scalar tij = this->transmissibility_(domainI,
-                                                           domainJ,
-                                                           insideDistance,
-                                                           outsideDistance,
-                                                           diffusionCoefficient_(volVarsI, stokesPhaseIdx, compIdx) * volVarsI.extrusionFactor(),
-                                                           diffusionCoefficient_(volVarsJ, stokesPhaseIdx, compIdx) * volVarsJ.extrusionFactor(),
-                                                           diffCoeffAvgType);
-                diffusiveFlux[compIdx] += -avgMolarDensity * tij * deltaMoleFrac;
-            }
+            const int domainICompIdx = couplingCompIdx(domainI, compIdx);
+            const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
+
+            assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx));
+
+            const Scalar deltaMoleFrac = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIdx) - volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIdx);
+            const Scalar tij = this->transmissibility_(domainI,
+                                                       domainJ,
+                                                       insideDistance,
+                                                       outsideDistance,
+                                                       diffusionCoefficient_(volVarsI, couplingPhaseIdx(domainI), domainICompIdx) * volVarsI.extrusionFactor(),
+                                                       diffusionCoefficient_(volVarsJ, couplingPhaseIdx(domainJ), domainJCompIdx) * volVarsJ.extrusionFactor(),
+                                                       diffCoeffAvgType);
+            diffusiveFlux[domainICompIdx] += -avgMolarDensity * tij * deltaMoleFrac;
         }
 
         const Scalar cumulativeFlux = std::accumulate(diffusiveFlux.begin(), diffusiveFlux.end(), 0.0);
-        diffusiveFlux[stokesPhaseIdx] = -cumulativeFlux;
+        diffusiveFlux[couplingCompIdx(domainI, 0)] = -cumulativeFlux;
+
 
         if(!useMoles)
         {
             //convert everything to a mass flux
             for(int compIdx = 0; compIdx < numComponents; ++compIdx)
-                diffusiveFlux[compIdx] *= FluidSystem<darcyIdx>::molarMass(compIdx);
+                diffusiveFlux[compIdx] *= FluidSystem<i>::molarMass(compIdx);
         }
 
         return diffusiveFlux;
+        // return NumEqVector(0);
     }
 
     /*!
@@ -725,8 +821,8 @@ protected:
         const auto& outsideScv = (*scvs(outsideFvGeometry).begin());
 
         // convective fluxes
-        const Scalar insideTerm = insideVolVars.density(stokesPhaseIdx) * insideVolVars.enthalpy(stokesPhaseIdx);
-        const Scalar outsideTerm = outsideVolVars.density(stokesPhaseIdx) * outsideVolVars.enthalpy(stokesPhaseIdx);
+        const Scalar insideTerm = insideVolVars.density(couplingPhaseIdx(domainI)) * insideVolVars.enthalpy(couplingPhaseIdx(domainI));
+        const Scalar outsideTerm = outsideVolVars.density(couplingPhaseIdx(domainJ)) * outsideVolVars.enthalpy(couplingPhaseIdx(domainJ));
 
         flux += this->advectiveFlux(insideTerm, outsideTerm, velocity, insideIsUpstream);
 
@@ -736,16 +832,19 @@ protected:
 
         for (int compIdx = 0; compIdx < diffusiveFlux.size(); ++compIdx)
         {
-            const bool insideDiffFluxIsUpstream = std::signbit(diffusiveFlux[compIdx]);
+            const int domainICompIdx = couplingCompIdx(domainI, compIdx);
+            const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
+
+            const bool insideDiffFluxIsUpstream = std::signbit(diffusiveFlux[domainICompIdx]);
             const Scalar componentEnthalpy = insideDiffFluxIsUpstream ?
-                                             FluidSystem<i>::componentEnthalpy(insideVolVars.fluidState(), stokesPhaseIdx, compIdx)
-                                           : FluidSystem<j>::componentEnthalpy(outsideVolVars.fluidState(), stokesPhaseIdx, compIdx);
+                                             getComponentEnthalpy(insideVolVars, couplingPhaseIdx(domainI), domainICompIdx)
+                                           : getComponentEnthalpy(outsideVolVars, couplingPhaseIdx(domainJ), domainJCompIdx);
 
             // always use a mass-based calculation for the energy balance
             if (useMoles)
-                diffusiveFlux[compIdx] *= FluidSystem<i>::molarMass(compIdx);
+                diffusiveFlux[domainICompIdx] *= FluidSystem<i>::molarMass(domainICompIdx);
 
-            flux += diffusiveFlux[compIdx] * componentEnthalpy;
+            flux += diffusiveFlux[domainICompIdx] * componentEnthalpy;
         }
 
         return flux;
diff --git a/test/freeflow/ransnc/flatplatetestproblem.hh b/test/freeflow/ransnc/flatplatetestproblem.hh
index 040bb514dbbd222e7fbcf9e4ba123a4dbf7c830a..723868298f52de7cab857c6f6a591eaf9286d0cb 100644
--- a/test/freeflow/ransnc/flatplatetestproblem.hh
+++ b/test/freeflow/ransnc/flatplatetestproblem.hh
@@ -27,6 +27,7 @@
 #include <dune/grid/yaspgrid.hh>
 
 #include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 #include <dumux/freeflow/turbulenceproperties.hh>
 
@@ -81,16 +82,16 @@ namespace Properties
   #endif
 #endif
 
-NEW_PROP_TAG(FluidSystem);
-
-// Select the fluid system
-SET_TYPE_PROP(FlatPlateNCTestTypeTag, FluidSystem,
-              FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-SET_INT_PROP(FlatPlateNCTestTypeTag, PhaseIdx,
-             GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+// The fluid system
+SET_PROP(FlatPlateNCTestTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
-SET_INT_PROP(FlatPlateNCTestTypeTag, ReplaceCompEqIdx, GET_PROP_VALUE(TypeTag, PhaseIdx));
+// replace the main component balance eq with a total balance eq
+SET_INT_PROP(FlatPlateNCTestTypeTag, ReplaceCompEqIdx, 0);
 
 // Set the grid type
 SET_TYPE_PROP(FlatPlateNCTestTypeTag, Grid,
@@ -158,9 +159,8 @@ class FlatPlateNCTestProblem : public ZeroEqProblem<TypeTag>
     using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
 
     static constexpr auto dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
-    static const unsigned int phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr auto transportEqIdx = Indices::conti0EqIdx;
-    static constexpr auto transportCompIdx = Indices::conti0EqIdx;
+    static constexpr auto transportEqIdx = Indices::conti0EqIdx + 1;
+    static constexpr auto transportCompIdx = Indices::conti0EqIdx + 1;
 
 public:
     FlatPlateNCTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
@@ -171,6 +171,7 @@ public:
         FluidSystem::init();
         Dumux::TurbulenceProperties<Scalar, dimWorld, true> turbulenceProperties;
         FluidState fluidState;
+        const auto phaseIdx = 0;
         fluidState.setPressure(phaseIdx, 1e5);
         fluidState.setTemperature(temperature());
         fluidState.setMassFraction(phaseIdx, phaseIdx, 1.0);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh
index 2f16c830fa18853d5510ee20e21b4aa948ab8afb..f4d864dc0dcbda86f0408d097c90b6737d5869de 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/darcyproblem.hh
@@ -33,6 +33,7 @@
 
 #include "./../1pspatialparams.hh"
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 #include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
 
@@ -49,10 +50,12 @@ NEW_TYPE_TAG(DarcyOnePTwoCTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
 SET_TYPE_PROP(DarcyOnePTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
 
 // The fluid system
-SET_TYPE_PROP(DarcyOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// Use water as phase
-SET_INT_PROP(DarcyOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::liquidPhaseIdx);
+SET_PROP(DarcyOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Use moles
 SET_BOOL_PROP(DarcyOnePTwoCTypeTag, UseMoles, true);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh
index ebc528deaeec1b3c580a27e59f4888bd055fad7e..a3627d489667487e5fc4963024fe83218da11051 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/horizontalflow/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 
 #include <dumux/freeflow/navierstokes/problem.hh>
@@ -41,8 +42,13 @@ namespace Properties
 {
 NEW_TYPE_TAG(StokesOnePTwoCTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
 
-// Set the fluid system
-SET_TYPE_PROP(StokesOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+// The fluid system
+SET_PROP(StokesOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Set the grid type
 SET_TYPE_PROP(StokesOnePTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh
index c07e21fac319d11cd0cad8dc4fcccd34ab444781..6fc7a42f49a39447a355d41104ba10f2cae61332 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/darcyproblem.hh
@@ -33,6 +33,7 @@
 
 #include "./../1pspatialparams.hh"
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 #include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
 
@@ -49,10 +50,12 @@ NEW_TYPE_TAG(DarcyOnePTwoCTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
 SET_TYPE_PROP(DarcyOnePTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
 
 // The fluid system
-SET_TYPE_PROP(DarcyOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// Use water as phase
-SET_INT_PROP(DarcyOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::liquidPhaseIdx);
+SET_PROP(DarcyOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Use moles
 SET_BOOL_PROP(DarcyOnePTwoCTypeTag, UseMoles, true);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh
index 5d6fed0b9543207c88beb4ac1fa5763f23409cde..ade5e9b2e1f9d3c4168743c4297f0269d584eaf1 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_1p2c/verticalflow/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 
 #include <dumux/freeflow/navierstokes/problem.hh>
@@ -41,11 +42,13 @@ namespace Properties
 {
 NEW_TYPE_TAG(StokesOnePTwoCTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
 
-// Set the fluid system
-SET_TYPE_PROP(StokesOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// Use water as phase
-SET_INT_PROP(StokesOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::liquidPhaseIdx);
+// The fluid system
+SET_PROP(StokesOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::liquidPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 // Set the grid type
 SET_TYPE_PROP(StokesOnePTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/CMakeLists.txt b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/CMakeLists.txt
index 845360ae428b1349030650f3d696c2482b618545..d85a5533fb9ed362e8045b3f7d568aa2261d9ad6 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/CMakeLists.txt
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/CMakeLists.txt
@@ -6,9 +6,9 @@ dune_add_test(NAME test_stokes1p2cdarcy2p2chorizontal
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS       --script fuzzy
                              --files ${CMAKE_SOURCE_DIR}/test/references/test_stokes1p2cdarcy2p2chorizontal_stokes-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy2p2chorizontal_stokes-00039.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy2p2chorizontal_stokes-00040.vtu
                                      ${CMAKE_SOURCE_DIR}/test/references/test_stokes1p2cdarcy2p2chorizontal_darcy-reference.vtu
-                                     ${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy2p2chorizontal_darcy-00039.vtu
+                                     ${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy2p2chorizontal_darcy-00040.vtu
 
                              --command "${CMAKE_CURRENT_BINARY_DIR}/test_stokes1p2cdarcy2p2chorizontal test_stokes1p2cdarcy2p2chorizontal.input")
 
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh
index 7e0f64317079d5ebc22fd7973df772686dde3fe8..dd402578eb847da3a6a81bca50d4eb1c9d68ff88 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/h2oair.hh>
 
 #include <dumux/freeflow/navierstokes/problem.hh>
@@ -49,11 +50,13 @@ NEW_TYPE_TAG(StokesOnePTwoCTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, Navier
 // Set the grid type
 SET_TYPE_PROP(StokesOnePTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
 
-// set the fluid system
-SET_TYPE_PROP(StokesOnePTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// set phase index (air)
-SET_INT_PROP(StokesOnePTwoCTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+// The fluid system
+SET_PROP(StokesOnePTwoCTypeTag, FluidSystem)
+{
+  using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+  static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the water phase
+  using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
 
 SET_INT_PROP(StokesOnePTwoCTypeTag, ReplaceCompEqIdx, 3);
 
@@ -107,9 +110,6 @@ class StokesSubProblem : public NavierStokesProblem<TypeTag>
 
     static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
 
-    static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
-    static constexpr auto transportCompIdx = 1 - phaseIdx;
-
 public:
     StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
     : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
@@ -128,10 +128,6 @@ public:
      */
     // \{
 
-
-    bool shouldWriteRestartFile() const
-    { return false; }
-
    /*!
      * \brief Return the temperature within the domain in [K].
      */
@@ -178,18 +174,10 @@ public:
             values.setNeumann(Indices::conti0EqIdx + 1);
         }
 
-        if (onUpperBoundary_(globalPos))
-        {
-            values.setDirichlet(Indices::velocityXIdx);
-            values.setDirichlet(Indices::velocityYIdx);
-            values.setNeumann(Indices::conti0EqIdx);
-            values.setNeumann(Indices::conti0EqIdx + 1);
-        }
-
         if (onRightBoundary_(globalPos))
         {
             values.setDirichlet(Indices::pressureIdx);
-            values.setOutflow(Indices::conti0EqIdx);
+            values.setOutflow(Indices::conti0EqIdx + 1);
 
 #if NONISOTHERMAL
             values.setOutflow(Indices::energyBalanceIdx);
@@ -242,17 +230,17 @@ public:
         FluidState fluidState;
         updateFluidStateForBC_(fluidState, elemVolVars[scv].pressure());
 
-        const Scalar density = useMoles ? fluidState.molarDensity(phaseIdx) : fluidState.density(phaseIdx);
+        const Scalar density = useMoles ? fluidState.molarDensity(0) : fluidState.density(0);
         const Scalar xVelocity = xVelocity_(globalPos);
 
         if (onLeftBoundary_(globalPos))
         {
             // rho*v*X at inflow
-            values[Indices::conti0EqIdx] = -xVelocity * density * refMoleFrac();
-            values[Indices::conti0EqIdx + 1] = -xVelocity * density * (1.0 - refMoleFrac());
+            values[Indices::conti0EqIdx + 1] = -xVelocity * density * refMoleFrac();
+            values[Indices::conti0EqIdx] = -xVelocity * density * (1.0 - refMoleFrac());
 
 #if NONISOTHERMAL
-            values[Indices::energyBalanceIdx] = -xVelocity * fluidState.density(phaseIdx) * fluidState.enthalpy(phaseIdx);
+            values[Indices::energyBalanceIdx] = -xVelocity * fluidState.density(0) * fluidState.enthalpy(0);
 #endif
         }
 
@@ -301,11 +289,11 @@ public:
         FluidState fluidState;
         updateFluidStateForBC_(fluidState, refPressure());
 
-        const Scalar density = FluidSystem::density(fluidState, phaseIdx);
+        const Scalar density = FluidSystem::density(fluidState, 0);
 
         PrimaryVariables values(0.0);
         values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
-        values[Indices::conti0EqIdx] = refMoleFrac();
+        values[Indices::conti0EqIdx + 1] = refMoleFrac();
         values[Indices::velocityXIdx] = xVelocity_(globalPos);
 
 #if NONISOTHERMAL
@@ -373,22 +361,22 @@ private:
     void updateFluidStateForBC_(FluidState& fluidState, const Scalar pressure) const
     {
         fluidState.setTemperature(refTemperature());
-        fluidState.setPressure(phaseIdx, pressure);
-        fluidState.setSaturation(phaseIdx, 1.0);
-        fluidState.setMoleFraction(phaseIdx, transportCompIdx, refMoleFrac());
-        fluidState.setMoleFraction(phaseIdx, phaseIdx, 1.0 - refMoleFrac());
+        fluidState.setPressure(0, pressure);
+        fluidState.setSaturation(0, 1.0);
+        fluidState.setMoleFraction(0, 1, refMoleFrac());
+        fluidState.setMoleFraction(0, 0, 1.0 - refMoleFrac());
 
         typename FluidSystem::ParameterCache paramCache;
-        paramCache.updatePhase(fluidState, phaseIdx);
+        paramCache.updatePhase(fluidState, 0);
 
-        const Scalar density = FluidSystem::density(fluidState, paramCache, phaseIdx);
-        fluidState.setDensity(phaseIdx, density);
+        const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+        fluidState.setDensity(0, density);
 
-        const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
-        fluidState.setMolarDensity(phaseIdx, molarDensity);
+        const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+        fluidState.setMolarDensity(0, molarDensity);
 
-        const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
-        fluidState.setEnthalpy(phaseIdx, enthalpy);
+        const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+        fluidState.setEnthalpy(0, enthalpy);
     }
 
     //! \brief set the profile of the inflow velocity (horizontal direction)
diff --git a/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh b/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh
index df3a7290e5e3c24fce2f646d6706fe4b515785e2..f4264647660c22e656d054b04be3066b15125d62 100644
--- a/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p_2p/stokesproblem.hh
@@ -26,6 +26,7 @@
 
 #include <dune/grid/yaspgrid.hh>
 
+#include <dumux/material/fluidsystems/1padapter.hh>
 #include <dumux/material/fluidsystems/2pimmiscible.hh>
 #include <dumux/material/fluidsystems/1pliquid.hh>
 #include <dumux/material/fluidsystems/1pgas.hh>
@@ -53,16 +54,16 @@ SET_TYPE_PROP(StokesOnePTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffset
 SET_PROP(StokesOnePTypeTag, FluidSystem)
 {
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using H2OType = Dumux::Components::SimpleH2O<Scalar>;
-    using H2OPhase = Dumux::FluidSystems::OnePLiquid<Scalar, H2OType>;
-    using AirType = Dumux::Components::TabulatedComponent<Components::Air<Scalar>, false >;
-    using AirPhase = Dumux::FluidSystems::OnePGas<Scalar, AirType>;
-    using type = FluidSystems::TwoPImmiscible<Scalar, H2OPhase, AirPhase> ;
+    using H2OType = Components::SimpleH2O<Scalar>;
+    using H2OPhase = FluidSystems::OnePLiquid<Scalar, H2OType>;
+    using AirType = Components::TabulatedComponent<Components::Air<Scalar>, false >;
+    using AirPhase = FluidSystems::OnePGas<Scalar, AirType>;
+    using TwoPFluidSystem = FluidSystems::TwoPImmiscible<Scalar, H2OPhase, AirPhase> ;
+
+    static constexpr auto phaseIdx = 1; // simulate the air phase
+    using type = FluidSystems::OnePAdapter<TwoPFluidSystem, phaseIdx>;
 };
 
-// consider air as phase
-SET_INT_PROP(StokesOnePTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::phase1Idx);
-
 // Set the problem property
 SET_TYPE_PROP(StokesOnePTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
 
diff --git a/test/references/test_stokes1p2cdarcy2p2chorizontal_darcy-reference.vtu b/test/references/test_stokes1p2cdarcy2p2chorizontal_darcy-reference.vtu
index 2f66edddb19d6e603e511f116a807be36fe58453..3ff49586547e583fbee4d8fa4aea958f0904f35a 100644
--- a/test/references/test_stokes1p2cdarcy2p2chorizontal_darcy-reference.vtu
+++ b/test/references/test_stokes1p2cdarcy2p2chorizontal_darcy-reference.vtu
@@ -4,18 +4,18 @@
     <Piece NumberOfCells="64" NumberOfPoints="81">
       <CellData Scalars="S_n" Vectors="velocity_liquid (m/s)">
         <DataArray type="Float32" Name="S_n" NumberOfComponents="1" format="ascii">
-          0.377422 0.377422 0.377422 0.377421 0.377421 0.377421 0.377421 0.377421 0.87225 0.87225 0.87225 0.87225
-          0.87225 0.87225 0.87225 0.87225 0.953604 0.953603 0.953602 0.953602 0.953601 0.953601 0.953601 0.9536
-          0.973905 0.973901 0.973896 0.973893 0.97389 0.973887 0.973885 0.973883 0.982539 0.982522 0.982509 0.982498
-          0.982489 0.982481 0.982475 0.982469 0.98896 0.988901 0.988857 0.988822 0.988793 0.988768 0.988747 0.988728
+          0.385645 0.385645 0.385645 0.385645 0.385645 0.385644 0.385644 0.385644 0.874632 0.874632 0.874632 0.874632
+          0.874632 0.874631 0.874631 0.874631 0.954314 0.954313 0.954312 0.954312 0.954311 0.954311 0.954311 0.95431
+          0.974243 0.974239 0.974234 0.97423 0.974227 0.974225 0.974223 0.974221 0.982745 0.982728 0.982714 0.982703
+          0.982694 0.982686 0.98268 0.982674 0.989175 0.989114 0.989068 0.989032 0.989002 0.988976 0.988954 0.988935
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="S_w" NumberOfComponents="1" format="ascii">
-          0.622578 0.622578 0.622578 0.622579 0.622579 0.622579 0.622579 0.622579 0.12775 0.12775 0.12775 0.12775
-          0.12775 0.12775 0.12775 0.12775 0.0463964 0.046397 0.0463977 0.0463982 0.0463986 0.046399 0.0463993 0.0463995
-          0.0260946 0.0260995 0.0261039 0.0261075 0.0261104 0.0261129 0.0261151 0.0261168 0.0174613 0.0174779 0.0174915 0.0175022
-          0.0175111 0.0175188 0.0175255 0.0175308 0.0110396 0.0110987 0.0111434 0.0111782 0.011207 0.0112318 0.0112535 0.0112718
+          0.614355 0.614355 0.614355 0.614355 0.614355 0.614356 0.614356 0.614356 0.125368 0.125368 0.125368 0.125368
+          0.125368 0.125369 0.125369 0.125369 0.0456864 0.045687 0.0456877 0.0456882 0.0456887 0.045689 0.0456894 0.0456896
+          0.0257565 0.0257615 0.025766 0.0257696 0.0257726 0.0257752 0.0257774 0.0257791 0.0172554 0.0172722 0.017286 0.0172969
+          0.0173059 0.0173137 0.0173205 0.0173259 0.0108245 0.0108859 0.0109323 0.0109684 0.0109982 0.0110238 0.0110463 0.0110652
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0
         </DataArray>
@@ -28,18 +28,18 @@
           100000 100000 100000 100000
         </DataArray>
         <DataArray type="Float32" Name="p_w" NumberOfComponents="1" format="ascii">
-          98503.8 98503.8 98503.8 98503.8 98503.8 98503.8 98503.8 98503.8 97797 97797 97797 97797
-          97797 97797 97797 97797 97324.7 97324.7 97324.7 97324.7 97324.7 97324.7 97324.7 97324.7
-          96994.4 96994.6 96994.7 96994.8 96994.8 96994.9 96994.9 96995 96711.2 96711.9 96712.6 96713
-          96713.4 96713.8 96714.1 96714.3 96351.1 96354.6 96357.2 96359.3 96361 96362.4 96363.7 96364.8
+          98495.8 98495.8 98495.8 98495.8 98495.8 98495.8 98495.8 98495.8 97788.9 97788.9 97789 97789
+          97789 97789 97789 97789 97316.6 97316.7 97316.7 97316.7 97316.7 97316.7 97316.7 97316.7
+          96986.1 96986.2 96986.3 96986.4 96986.5 96986.6 96986.6 96986.7 96701.8 96702.6 96703.2 96703.7
+          96704.1 96704.5 96704.8 96705.1 96338.3 96342 96344.7 96346.8 96348.6 96350.1 96351.5 96352.6
           95697.6 95697.6 95697.6 95697.6 95697.6 95697.6 95697.6 95697.6 95697.5 95697.5 95697.5 95697.5
           95697.5 95697.5 95697.5 95697.5
         </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
-          1498.52 1498.52 1498.52 1498.52 1498.52 1498.52 1498.51 1498.51 2204.53 2204.53 2204.53 2204.53
-          2204.53 2204.53 2204.53 2204.53 2676.26 2676.25 2676.24 2676.24 2676.23 2676.23 2676.22 2676.22
-          3006.16 3006.04 3005.93 3005.85 3005.77 3005.71 3005.66 3005.62 3289.16 3288.42 3287.8 3287.32
-          3286.92 3286.58 3286.28 3286.04 3649.14 3645.64 3642.99 3640.93 3639.23 3637.76 3636.48 3635.4
+          1506.54 1506.54 1506.54 1506.54 1506.54 1506.54 1506.54 1506.54 2212.56 2212.56 2212.56 2212.56
+          2212.56 2212.56 2212.56 2212.56 2684.31 2684.31 2684.3 2684.29 2684.29 2684.28 2684.28 2684.28
+          3014.5 3014.37 3014.26 3014.17 3014.09 3014.03 3013.98 3013.93 3298.52 3297.75 3297.11 3296.62
+          3296.2 3295.85 3295.54 3295.29 3661.86 3658.23 3655.49 3653.35 3651.58 3650.07 3648.74 3647.62
           4302.49 4302.49 4302.49 4302.49 4302.49 4302.49 4302.49 4302.49 4302.49 4302.49 4302.49 4302.49
           4302.49 4302.49 4302.49 4302.49
         </DataArray>
@@ -48,96 +48,96 @@
           997.051 997.051 997.051 997.051 997.051 997.051 997.051 997.051 997.051 997.051 997.051 997.051
           997.051 997.051 997.051 997.051 997.051 997.051 997.051 997.051 997.05 997.05 997.05 997.05
           997.05 997.05 997.05 997.05 997.05 997.05 997.05 997.05 997.05 997.05 997.05 997.05
-          628.517 647.342 661.237 672.15 681.269 689.17 696.17 702.234 140.491 165.19 182.192 195.534
-          206.696 216.37 224.954 232.696
+          608.982 627.135 640.529 651.046 659.835 667.451 674.197 680.04 136.068 159.972 176.422 189.327
+          200.123 209.48 217.781 225.261
         </DataArray>
         <DataArray type="Float32" Name="rho_n" NumberOfComponents="1" format="ascii">
           1.15428 1.15428 1.15428 1.15428 1.15428 1.15428 1.15428 1.15428 1.15427 1.15427 1.15427 1.15427
           1.15427 1.15427 1.15427 1.15427 1.15427 1.15427 1.15427 1.15427 1.15427 1.15427 1.15427 1.15427
           1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426
           1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426 1.15426
-          1.15941 1.15915 1.15896 1.1588 1.15868 1.15856 1.15847 1.15838 1.16626 1.16591 1.16567 1.16549
-          1.16533 1.16519 1.16507 1.16496
+          1.15969 1.15943 1.15925 1.1591 1.15898 1.15887 1.15877 1.15869 1.16632 1.16598 1.16575 1.16557
+          1.16542 1.16529 1.16517 1.16507
         </DataArray>
         <DataArray type="Float32" Name="mob_w" NumberOfComponents="1" format="ascii">
-          243.306 243.306 243.306 243.306 243.306 243.307 243.307 243.307 2.25177 2.25177 2.25178 2.25178
-          2.25179 2.25179 2.2518 2.2518 0.101914 0.101919 0.101923 0.101927 0.10193 0.101932 0.101935 0.101936
-          0.0150021 0.0150119 0.0150209 0.0150281 0.0150341 0.0150392 0.0150436 0.0150469 0.00336369 0.00337637 0.00338686 0.00339514
-          0.00340202 0.00340794 0.00341311 0.00341725 0.000430107 0.000442158 0.000451431 0.000458737 0.000464837 0.000470119 0.000474785 0.000478729
+          233.688 233.688 233.688 233.688 233.688 233.689 233.689 233.689 2.1293 2.1293 2.1293 2.12931
+          2.12931 2.12932 2.12932 2.12932 0.0970225 0.097027 0.0970313 0.0970349 0.097038 0.0970405 0.0970427 0.0970443
+          0.0143291 0.0143388 0.0143477 0.0143548 0.0143607 0.0143657 0.01437 0.0143733 0.00320822 0.00322071 0.00323103 0.00323917
+          0.00324593 0.00325174 0.00325682 0.00326089 0.000388031 0.000399745 0.000408758 0.000415861 0.000421793 0.000426928 0.000431465 0.000435299
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0
         </DataArray>
         <DataArray type="Float32" Name="mob_n" NumberOfComponents="1" format="ascii">
-          9036.5 9036.5 9036.49 9036.48 9036.47 9036.46 9036.46 9036.45 45130.6 45130.6 45130.5 45130.5
-          45130.5 45130.5 45130.5 45130.5 53066.8 53066.8 53066.7 53066.7 53066.6 53066.6 53066.6 53066.5
-          54566.3 54566 54565.8 54565.6 54565.4 54565.3 54565.1 54565.1 54956.9 54956.4 54955.9 54955.5
-          54955.2 54954.9 54954.7 54954.5 55127 55125.9 55125.1 55124.5 55123.9 55123.5 55123.1 55122.8
-          54850.3 54867.2 54879.6 54889.4 54897.5 54904.6 54910.9 54916.3 54417.2 54438.9 54453.9 54465.7
-          54475.6 54484.1 54491.7 54498.6
+          9439.38 9439.38 9439.37 9439.36 9439.35 9439.34 9439.34 9439.33 45340.1 45340.1 45340.1 45340.1
+          45340.1 45340.1 45340.1 45340.1 53130.1 53130.1 53130 53130 53129.9 53129.9 53129.9 53129.9
+          54584.8 54584.5 54584.2 54584.1 54583.9 54583.8 54583.6 54583.5 54964.1 54963.5 54963 54962.6
+          54962.3 54962.1 54961.8 54961.6 55130.7 55129.7 55128.8 55128.2 55127.7 55127.2 55126.8 55126.5
+          54832.9 54849.1 54861.1 54870.5 54878.3 54885.2 54891.2 54896.4 54413.3 54434.3 54448.8 54460.2
+          54469.8 54478 54485.4 54492
         </DataArray>
         <DataArray type="Float32" Name="X^H2O_liquid" NumberOfComponents="1" format="ascii">
           0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
           0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
           0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
           0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978 0.999978
-          0.630352 0.649232 0.663169 0.674114 0.68326 0.691185 0.698205 0.704288 0.140878 0.16565 0.182702 0.196083
-          0.207278 0.216981 0.225591 0.233356
+          0.610759 0.628966 0.642399 0.652948 0.661763 0.669401 0.676167 0.682027 0.136441 0.160416 0.176914 0.189858
+          0.200687 0.210071 0.218396 0.225899
         </DataArray>
         <DataArray type="Float32" Name="X^Air_liquid" NumberOfComponents="1" format="ascii">
           2.15464e-05 2.15464e-05 2.15464e-05 2.15464e-05 2.15464e-05 2.15464e-05 2.15464e-05 2.15464e-05 2.15462e-05 2.15462e-05 2.15462e-05 2.15462e-05
           2.15462e-05 2.15462e-05 2.15462e-05 2.15462e-05 2.15461e-05 2.15461e-05 2.15461e-05 2.15461e-05 2.15461e-05 2.15461e-05 2.15461e-05 2.15461e-05
           2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05
           2.1546e-05 2.1546e-05 2.1546e-05 2.1546e-05 2.15459e-05 2.15459e-05 2.15459e-05 2.15459e-05 2.15459e-05 2.15459e-05 2.15459e-05 2.15459e-05
-          2.18069e-05 2.17935e-05 2.17837e-05 2.1776e-05 2.17695e-05 2.17639e-05 2.1759e-05 2.17547e-05 2.21515e-05 2.21342e-05 2.21223e-05 2.21129e-05
-          2.21051e-05 2.20983e-05 2.20922e-05 2.20867e-05
+          2.18207e-05 2.18079e-05 2.17984e-05 2.17909e-05 2.17847e-05 2.17793e-05 2.17745e-05 2.17704e-05 2.21546e-05 2.21379e-05 2.21264e-05 2.21173e-05
+          2.21097e-05 2.21031e-05 2.20973e-05 2.2092e-05
         </DataArray>
         <DataArray type="Float32" Name="X^H2O_gas" NumberOfComponents="1" format="ascii">
           0.0199835 0.0199835 0.0199835 0.0199835 0.0199835 0.0199835 0.0199835 0.0199835 0.0199837 0.0199837 0.0199837 0.0199837
           0.0199837 0.0199837 0.0199837 0.0199837 0.0199838 0.0199838 0.0199838 0.0199838 0.0199838 0.0199838 0.0199838 0.0199838
           0.0199839 0.0199839 0.0199839 0.0199839 0.0199839 0.0199839 0.0199839 0.0199839 0.0199839 0.0199839 0.0199839 0.0199839
           0.0199839 0.0199839 0.0199839 0.0199839 0.019984 0.019984 0.019984 0.019984 0.019984 0.019984 0.019984 0.019984
-          0.012541 0.0129196 0.0131991 0.0134188 0.0136023 0.0137614 0.0139024 0.0140245 0.00278644 0.00327737 0.00361547 0.00388087
-          0.00410299 0.00429555 0.00446644 0.0046206
+          0.0121483 0.0125132 0.0127825 0.0129941 0.0131709 0.0133242 0.01346 0.0135776 0.00269856 0.00317362 0.00350069 0.00375739
+          0.0039722 0.0041584 0.00432363 0.00447257
         </DataArray>
         <DataArray type="Float32" Name="X^Air_gas" NumberOfComponents="1" format="ascii">
           0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016
           0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016
           0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016
           0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016 0.980016
-          0.987459 0.98708 0.986801 0.986581 0.986398 0.986239 0.986098 0.985976 0.997214 0.996723 0.996385 0.996119
-          0.995897 0.995704 0.995534 0.995379
+          0.987852 0.987487 0.987217 0.987006 0.986829 0.986676 0.98654 0.986422 0.997301 0.996826 0.996499 0.996243
+          0.996028 0.995842 0.995676 0.995527
         </DataArray>
         <DataArray type="Float32" Name="x^H2O_liquid" NumberOfComponents="1" format="ascii">
           0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
           0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
           0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
           0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987 0.999987
-          0.63036 0.64924 0.663177 0.674122 0.683268 0.691193 0.698214 0.704296 0.140886 0.165658 0.182711 0.196092
-          0.207287 0.21699 0.225599 0.233364
+          0.610767 0.628974 0.642407 0.652956 0.661771 0.669409 0.676175 0.682036 0.13645 0.160425 0.176923 0.189867
+          0.200695 0.210079 0.218405 0.225908
         </DataArray>
         <DataArray type="Float32" Name="x^Air_liquid" NumberOfComponents="1" format="ascii">
           1.34035e-05 1.34035e-05 1.34035e-05 1.34035e-05 1.34035e-05 1.34035e-05 1.34035e-05 1.34035e-05 1.34034e-05 1.34034e-05 1.34034e-05 1.34034e-05
           1.34034e-05 1.34034e-05 1.34034e-05 1.34034e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05
           1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34033e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05
           1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05 1.34032e-05
-          1.35656e-05 1.35573e-05 1.35512e-05 1.35464e-05 1.35424e-05 1.35389e-05 1.35358e-05 1.35331e-05 1.37807e-05 1.37698e-05 1.37623e-05 1.37564e-05
-          1.37515e-05 1.37472e-05 1.37434e-05 1.374e-05
+          1.35742e-05 1.35662e-05 1.35603e-05 1.35557e-05 1.35518e-05 1.35484e-05 1.35455e-05 1.35429e-05 1.37826e-05 1.37721e-05 1.37648e-05 1.37591e-05
+          1.37544e-05 1.37502e-05 1.37466e-05 1.37433e-05
         </DataArray>
         <DataArray type="Float32" Name="x^H2O_gas" NumberOfComponents="1" format="ascii">
           0.0317389 0.0317389 0.0317389 0.0317389 0.0317389 0.0317389 0.0317389 0.0317389 0.0317391 0.0317391 0.0317391 0.0317391
           0.0317391 0.0317391 0.0317391 0.0317391 0.0317393 0.0317393 0.0317393 0.0317393 0.0317393 0.0317393 0.0317393 0.0317393
           0.0317394 0.0317394 0.0317394 0.0317394 0.0317394 0.0317394 0.0317394 0.0317394 0.0317395 0.0317395 0.0317395 0.0317395
           0.0317395 0.0317395 0.0317395 0.0317395 0.0317396 0.0317396 0.0317396 0.0317396 0.0317396 0.0317396 0.0317396 0.0317396
-          0.0200076 0.0206069 0.0210493 0.0213967 0.021687 0.0219385 0.0221613 0.0223544 0.00447173 0.005258 0.00579925 0.00622395
-          0.00657929 0.00688726 0.00716052 0.00740698
+          0.0193857 0.0199637 0.02039 0.0207248 0.0210046 0.0212471 0.0214618 0.0216478 0.00433092 0.00509189 0.00561554 0.00602637
+          0.00637007 0.00666792 0.00693217 0.00717032
         </DataArray>
         <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
           0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261
           0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261
           0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.968261 0.96826 0.96826 0.96826 0.96826
           0.96826 0.96826 0.96826 0.96826 0.96826 0.96826 0.96826 0.96826 0.96826 0.96826 0.96826 0.96826
-          0.979992 0.979393 0.978951 0.978603 0.978313 0.978061 0.977839 0.977646 0.995528 0.994742 0.994201 0.993776
-          0.993421 0.993113 0.992839 0.992593
+          0.980614 0.980036 0.97961 0.979275 0.978995 0.978753 0.978538 0.978352 0.995669 0.994908 0.994384 0.993974
+          0.99363 0.993332 0.993068 0.99283
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41 0.41
@@ -156,40 +156,40 @@
           2 2 2 2
         </DataArray>
         <DataArray type="Float32" Name="velocity_liquid (m/s)" NumberOfComponents="3" format="ascii">
-          -8.48001e-11 1.09813e-08 0 -1.12443e-10 1.08994e-08 0 -1.52919e-10 1.08172e-08 0 -1.66394e-10 1.07496e-08 0
-          -1.57571e-10 1.06935e-08 0 -1.30033e-10 1.06457e-08 0 -8.65701e-11 1.06051e-08 0 -6.11828e-11 1.05762e-08 0
-          -2.53782e-12 1.20537e-08 0 -2.79664e-12 1.19699e-08 0 -3.01908e-12 1.18872e-08 0 -2.86446e-12 1.18194e-08 0
-          -2.57002e-12 1.17632e-08 0 -2.14933e-12 1.17153e-08 0 -1.54464e-12 1.16745e-08 0 -1.18445e-12 1.16452e-08 0
-          -6.09198e-12 1.33422e-08 0 -6.02142e-12 1.32527e-08 0 -5.41581e-12 1.31688e-08 0 -4.46749e-12 1.3101e-08 0
-          -3.7586e-12 1.30448e-08 0 -3.20619e-12 1.29968e-08 0 -2.53963e-12 1.29556e-08 0 -2.12985e-12 1.29251e-08 0
-          -1.51165e-11 1.36316e-08 0 -1.44346e-11 1.35314e-08 0 -1.24007e-11 1.34464e-08 0 -1.00983e-11 1.33791e-08 0
-          -8.49103e-12 1.33233e-08 0 -7.29549e-12 1.32755e-08 0 -5.92978e-12 1.32338e-08 0 -5.10282e-12 1.3201e-08 0
-          -2.12609e-11 1.37378e-08 0 -1.94163e-11 1.3624e-08 0 -1.57246e-11 1.35384e-08 0 -1.2693e-11 1.34717e-08 0
-          -1.07042e-11 1.34164e-08 0 -9.27515e-12 1.33687e-08 0 -7.79459e-12 1.33267e-08 0 -6.93867e-12 1.32911e-08 0
-          -1.31054e-11 1.03443e-08 0 -1.16197e-11 1.0397e-08 0 -9.07607e-12 1.04416e-08 0 -7.36891e-12 1.04776e-08 0
-          -6.27869e-12 1.05079e-08 0 -5.505e-12 1.05343e-08 0 -4.77771e-12 1.05576e-08 0 -4.38302e-12 1.05767e-08 0
-          0 3.45894e-09 0 0 3.57197e-09 0 0 3.65944e-09 0 0 3.72858e-09 0
-          0 3.78644e-09 0 0 3.83664e-09 0 0 3.88105e-09 0 0 3.91861e-09 0
+          -8.35979e-11 1.06976e-08 0 -1.10834e-10 1.06168e-08 0 -1.50703e-10 1.05358e-08 0 -1.63955e-10 1.0469e-08 0
+          -1.55237e-10 1.04138e-08 0 -1.28088e-10 1.03668e-08 0 -8.52653e-11 1.03268e-08 0 -6.02559e-11 1.02984e-08 0
+          -2.46496e-12 1.17102e-08 0 -2.7154e-12 1.16275e-08 0 -2.92989e-12 1.15459e-08 0 -2.77879e-12 1.14791e-08 0
+          -2.49242e-12 1.14238e-08 0 -2.08385e-12 1.13766e-08 0 -1.49711e-12 1.13364e-08 0 -1.14772e-12 1.13076e-08 0
+          -6.04702e-12 1.29272e-08 0 -5.97586e-12 1.28389e-08 0 -5.37252e-12 1.27563e-08 0 -4.42943e-12 1.26894e-08 0
+          -3.72468e-12 1.26341e-08 0 -3.17573e-12 1.25869e-08 0 -2.51399e-12 1.25463e-08 0 -2.10739e-12 1.25163e-08 0
+          -1.49736e-11 1.32017e-08 0 -1.42933e-11 1.31028e-08 0 -1.22718e-11 1.3019e-08 0 -9.98795e-12 1.29526e-08 0
+          -8.39427e-12 1.28977e-08 0 -7.20914e-12 1.28506e-08 0 -5.85639e-12 1.28096e-08 0 -5.03765e-12 1.27774e-08 0
+          -2.09959e-11 1.33032e-08 0 -1.91659e-11 1.31908e-08 0 -1.55097e-11 1.31064e-08 0 -1.2513e-11 1.30407e-08 0
+          -1.05478e-11 1.29862e-08 0 -9.13602e-12 1.29393e-08 0 -7.67388e-12 1.28979e-08 0 -6.82881e-12 1.28629e-08 0
+          -1.23029e-11 9.72414e-09 0 -1.091e-11 9.77233e-09 0 -8.52514e-12 9.81355e-09 0 -6.92447e-12 9.84682e-09 0
+          -5.90201e-12 9.87489e-09 0 -5.17611e-12 9.89938e-09 0 -4.49226e-12 9.92105e-09 0 -4.12061e-12 9.93867e-09 0
+          0 3.05673e-09 0 0 3.16449e-09 0 0 3.24791e-09 0 0 3.31388e-09 0
+          0 3.36911e-09 0 0 3.41703e-09 0 0 3.45944e-09 0 0 3.49529e-09 0
           0 -0 0 0 -0 0 0 -0 0 0 -0 0
           0 -0 0 0 -0 0 0 -0 0 0 -0 0
         </DataArray>
         <DataArray type="Float32" Name="velocity_gas (m/s)" NumberOfComponents="3" format="ascii">
-          3.27518e-09 2.72687e-09 0 4.64402e-09 3.719e-09 0 6.91811e-09 5.36027e-09 0 8.1506e-09 7.23489e-09 0
-          8.17323e-09 9.00974e-09 0 6.97485e-09 1.04594e-08 0 4.70258e-09 1.14981e-08 0 3.32408e-09 1.20484e-08 0
-          2.7815e-08 -3.54016e-08 0 3.83997e-08 -2.79935e-08 0 5.48529e-08 -1.71205e-08 0 6.17797e-08 -5.61871e-09 0
-          5.94391e-08 5.21354e-09 0 4.92032e-08 1.39725e-08 0 3.26612e-08 2.11668e-08 0 2.29563e-08 2.5521e-08 0
-          5.40177e-08 -1.07033e-07 0 7.13046e-08 -8.12385e-08 0 9.62818e-08 -4.93371e-08 0 1.04453e-07 -1.93708e-08 0
-          9.99779e-08 7.28965e-09 0 8.45185e-08 3.01663e-08 0 5.78984e-08 5.0419e-08 0 4.17814e-08 6.44544e-08 0
-          9.43572e-08 -1.79382e-07 0 1.15181e-07 -1.19947e-07 0 1.41397e-07 -6.35624e-08 0 1.44714e-07 -1.86248e-08 0
-          1.35857e-07 1.86674e-08 0 1.16939e-07 5.27508e-08 0 8.44036e-08 8.66973e-08 0 6.40043e-08 1.16149e-07 0
-          1.50602e-07 -2.59527e-07 0 1.67452e-07 -1.46384e-07 0 1.83567e-07 -6.7674e-08 0 1.77402e-07 -1.40306e-08 0
-          1.64263e-07 2.83716e-08 0 1.45107e-07 6.88993e-08 0 1.12515e-07 1.14955e-07 0 9.13668e-08 1.67394e-07 0
-          2.17382e-07 4.38789e-06 0 2.21304e-07 4.32928e-06 0 2.17149e-07 4.24336e-06 0 2.00779e-07 4.15937e-06 0
-          1.84513e-07 4.0857e-06 0 1.67075e-07 4.02727e-06 0 1.39723e-07 3.99178e-06 0 1.21835e-07 3.99548e-06 0
-          1.31117e-07 1.38062e-05 0 1.36116e-07 1.37137e-05 0 1.38999e-07 1.35752e-05 0 1.33942e-07 1.3446e-05 0
-          1.28041e-07 1.33344e-05 0 1.2142e-07 1.3244e-05 0 1.09776e-07 1.31837e-05 0 1.01791e-07 1.31871e-05 0
-          1.35717e-07 1.85359e-05 0 1.32641e-07 1.86131e-05 0 1.26482e-07 1.85964e-05 0 1.21127e-07 1.85623e-05 0
-          1.17324e-07 1.853e-05 0 1.14703e-07 1.85087e-05 0 1.12329e-07 1.85103e-05 0 1.11045e-07 1.85757e-05 0
+          3.62987e-09 3.56115e-09 0 5.07787e-09 4.43266e-09 0 7.45539e-09 5.92228e-09 0 8.69744e-09 7.64949e-09 0
+          8.66543e-09 9.29752e-09 0 7.36466e-09 1.06501e-08 0 4.95193e-09 1.16246e-08 0 3.49543e-09 1.21437e-08 0
+          2.71153e-08 -3.46414e-08 0 3.75032e-08 -2.73448e-08 0 5.36754e-08 -1.6622e-08 0 6.05211e-08 -5.255e-09 0
+          5.82867e-08 5.4664e-09 0 4.82983e-08 1.41595e-08 0 3.20806e-08 2.12995e-08 0 2.25554e-08 2.56211e-08 0
+          5.29274e-08 -1.05649e-07 0 6.99057e-08 -8.02088e-08 0 9.44718e-08 -4.86921e-08 0 1.02586e-07 -1.9024e-08 0
+          9.82577e-08 7.40394e-09 0 8.30966e-08 3.01316e-08 0 5.69483e-08 5.02424e-08 0 4.11067e-08 6.41673e-08 0
+          9.25658e-08 -1.76561e-07 0 1.1305e-07 -1.1819e-07 0 1.3889e-07 -6.26746e-08 0 1.42264e-07 -1.83177e-08 0
+          1.33654e-07 1.8557e-08 0 1.15113e-07 5.23179e-08 0 8.31304e-08 8.594e-08 0 6.30593e-08 1.15086e-07 0
+          1.47737e-07 -2.55183e-07 0 1.64379e-07 -1.4423e-07 0 1.80394e-07 -6.68122e-08 0 1.74508e-07 -1.38885e-08 0
+          1.6172e-07 2.80414e-08 0 1.42966e-07 6.81678e-08 0 1.10922e-07 1.13753e-07 0 9.01036e-08 1.65603e-07 0
+          2.13192e-07 4.64662e-06 0 2.17223e-07 4.59281e-06 0 2.13449e-07 4.51169e-06 0 1.97583e-07 4.43194e-06 0
+          1.81745e-07 4.36196e-06 0 1.64698e-07 4.30672e-06 0 1.37822e-07 4.27386e-06 0 1.20215e-07 4.2794e-06 0
+          1.30872e-07 1.37668e-05 0 1.35744e-07 1.36795e-05 0 1.38477e-07 1.35465e-05 0 1.3341e-07 1.34222e-05 0
+          1.27548e-07 1.33148e-05 0 1.20993e-07 1.32281e-05 0 1.09461e-07 1.31708e-05 0 1.01549e-07 1.31763e-05 0
+          1.35479e-07 1.79447e-05 0 1.32417e-07 1.80199e-05 0 1.26279e-07 1.80031e-05 0 1.20939e-07 1.79695e-05 0
+          1.17147e-07 1.79379e-05 0 1.14535e-07 1.79171e-05 0 1.12177e-07 1.79189e-05 0 1.10905e-07 1.79841e-05 0
         </DataArray>
         <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
diff --git a/test/references/test_stokes1p2cdarcy2p2chorizontal_stokes-reference.vtu b/test/references/test_stokes1p2cdarcy2p2chorizontal_stokes-reference.vtu
index e49f9e0f6129c7b503ed597fdae6e5b28e881f4f..71f578ca3e976f25959c248947f73ca311b9c521 100644
--- a/test/references/test_stokes1p2cdarcy2p2chorizontal_stokes-reference.vtu
+++ b/test/references/test_stokes1p2cdarcy2p2chorizontal_stokes-reference.vtu
@@ -2,7 +2,7 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="128" NumberOfPoints="153">
-      <CellData Scalars="p" Vectors="velocity_liquid (m/s)">
+      <CellData Scalars="p" Vectors="velocity_gas (m/s)">
         <DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
           100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
           100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
@@ -17,8 +17,8 @@
           99997.6 99997.6 99997.6 99997.6 99997.6 99997.6 99997.6 99997.6
         </DataArray>
         <DataArray type="Float32" Name="rhoMolar" NumberOfComponents="1" format="ascii">
-          40.3395 40.3395 40.3395 40.3395 40.3396 40.3396 40.3396 40.3396 40.3395 40.3395 40.3395 40.3395
-          40.3395 40.3396 40.3396 40.3396 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3396
+          40.3395 40.3395 40.3395 40.3395 40.3395 40.3396 40.3396 40.3396 40.3395 40.3395 40.3395 40.3395
+          40.3395 40.3395 40.3396 40.3396 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3396
           40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395
           40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395
           40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395 40.3395
@@ -30,43 +30,69 @@
           40.3385 40.3385 40.3385 40.3385 40.3385 40.3385 40.3385 40.3385
         </DataArray>
         <DataArray type="Float32" Name="rho" NumberOfComponents="1" format="ascii">
-          1.16714 1.16678 1.16654 1.16634 1.16618 1.16605 1.16592 1.16581 1.16722 1.16686 1.16661 1.16642
-          1.16626 1.16612 1.166 1.16589 1.16734 1.16698 1.16673 1.16654 1.16638 1.16624 1.16611 1.166
-          1.1675 1.16715 1.1669 1.16671 1.16655 1.16641 1.16629 1.16617 1.16772 1.16739 1.16715 1.16696
-          1.1668 1.16666 1.16654 1.16642 1.16795 1.16769 1.16748 1.16731 1.16715 1.16702 1.1669 1.16678
-          1.16814 1.168 1.16786 1.16772 1.16759 1.16748 1.16737 1.16726 1.16822 1.16818 1.16813 1.16807
-          1.16801 1.16794 1.16787 1.16779 1.16823 1.16823 1.16822 1.16821 1.1682 1.16818 1.16816 1.16814
-          1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16822 1.16822 1.16823 1.16823 1.16823 1.16823
+          1.16717 1.16683 1.16659 1.1664 1.16625 1.16611 1.166 1.16589 1.16725 1.1669 1.16667 1.16648
+          1.16632 1.16619 1.16607 1.16596 1.16736 1.16702 1.16678 1.16659 1.16644 1.1663 1.16618 1.16607
+          1.16752 1.16718 1.16695 1.16676 1.1666 1.16647 1.16635 1.16624 1.16773 1.16742 1.16719 1.167
+          1.16685 1.16671 1.16659 1.16648 1.16796 1.16771 1.16751 1.16733 1.16719 1.16706 1.16694 1.16683
+          1.16814 1.16801 1.16787 1.16774 1.16761 1.1675 1.16739 1.16729 1.16822 1.16818 1.16814 1.16808
+          1.16801 1.16795 1.16788 1.16781 1.16823 1.16823 1.16822 1.16821 1.1682 1.16818 1.16817 1.16814
+          1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16822 1.16823 1.16823 1.16823 1.16823
           1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16823 1.16823
           1.16822 1.16822 1.16822 1.16822 1.16822 1.16822 1.16822 1.16822 1.16822 1.16822 1.16822 1.16822
           1.16822 1.16822 1.16822 1.16822 1.16821 1.16821 1.16821 1.16821 1.16821 1.16821 1.16821 1.16821
           1.1682 1.1682 1.1682 1.1682 1.1682 1.1682 1.1682 1.1682
         </DataArray>
+        <DataArray type="Float32" Name="X^Air_gas" NumberOfComponents="1" format="ascii">
+          0.998508 0.998017 0.997682 0.99742 0.997201 0.99701 0.996841 0.996688 0.998617 0.998124 0.99779 0.997527
+          0.997307 0.997116 0.996947 0.996793 0.998777 0.998284 0.997949 0.997686 0.997465 0.997274 0.997105 0.99695
+          0.999003 0.998519 0.998185 0.997922 0.997701 0.997509 0.997339 0.997184 0.999297 0.998848 0.998524 0.998264
+          0.998045 0.997854 0.997684 0.997529 0.999618 0.999267 0.998978 0.998735 0.998525 0.998341 0.998175 0.998022
+          0.999871 0.999685 0.99949 0.999302 0.999128 0.998967 0.998819 0.998679 0.99998 0.999934 0.999866 0.999784
+          0.999693 0.999598 0.999502 0.999403 0.999999 0.999995 0.999986 0.999973 0.999955 0.999934 0.999908 0.999878
+          1 1 0.999999 0.999999 0.999997 0.999996 0.999994 0.999991 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1
+        </DataArray>
+        <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
+          0.997604 0.996815 0.99628 0.995859 0.995507 0.995202 0.994932 0.994686 0.997779 0.996988 0.996451 0.99603
+          0.995678 0.995372 0.995101 0.994855 0.998035 0.997245 0.996707 0.996285 0.995932 0.995625 0.995354 0.995106
+          0.998398 0.997621 0.997086 0.996663 0.996309 0.996002 0.99573 0.995481 0.99887 0.99815 0.997629 0.997212
+          0.996861 0.996555 0.996283 0.996034 0.999387 0.998822 0.998358 0.997968 0.997631 0.997335 0.99707 0.996824
+          0.999793 0.999494 0.99918 0.998878 0.998598 0.998341 0.998103 0.997877 0.999968 0.999893 0.999785 0.999652
+          0.999506 0.999354 0.9992 0.999041 0.999998 0.999991 0.999978 0.999956 0.999928 0.999893 0.999852 0.999804
+          1 1 0.999999 0.999998 0.999996 0.999993 0.99999 0.999985 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1
+        </DataArray>
         <DataArray type="Float32" Name="X^H2O_gas" NumberOfComponents="1" format="ascii">
-          0.00154052 0.00204828 0.00239348 0.00266459 0.00289157 0.00308836 0.00326302 0.00342175 0.00142811 0.00193713 0.00228291 0.00255451
-          0.00278191 0.00297907 0.00315406 0.00331321 0.00126295 0.00177172 0.00211793 0.00239006 0.00261797 0.00281561 0.00299105 0.00315078
-          0.00102955 0.00152956 0.00187444 0.00214646 0.0023746 0.00257259 0.00274843 0.00290889 0.000726268 0.00118927 0.00152476 0.0017931
-          0.0020195 0.0022166 0.00239201 0.00255285 0.000394089 0.000756935 0.00105581 0.00130669 0.00152311 0.00171385 0.00188496 0.00204357
-          0.000133273 0.000324795 0.000527023 0.000720945 0.000900928 0.00106658 0.0012195 0.00136525 2.08739e-05 6.84652e-05 0.000138374 0.000223388
-          0.000317112 0.000414947 0.000514143 0.000616605 1.1804e-06 5.52265e-06 1.43412e-05 2.79917e-05 4.62355e-05 6.86364e-05 9.48696e-05 0.000126281
-          2.29444e-08 1.76069e-07 6.04067e-07 1.40896e-06 2.63735e-06 4.31097e-06 6.47027e-06 9.36497e-06 1.86331e-10 2.97042e-09 1.38611e-08 3.78557e-08
-          7.70538e-08 1.31579e-07 2.01878e-07 2.92495e-07 9.88026e-13 3.72966e-11 2.35447e-10 7.38633e-10 1.59414e-09 2.75281e-09 4.11998e-09 5.48111e-09
-          5.11072e-15 4.30405e-13 3.542e-12 1.25224e-11 2.83739e-11 4.92161e-11 7.11552e-11 8.64142e-11 0 5.11072e-15 5.1936e-14 2.01804e-13
-          4.75297e-13 8.2766e-13 1.16552e-12 1.34495e-12 0 0 8.28765e-16 3.31506e-15 8.0114e-15 1.38128e-14 1.91997e-14 2.1686e-14
-          0 0 -6.90638e-17 1.38128e-16 1.38128e-16 2.76255e-16 4.14383e-16 4.14383e-16
+          0.00149191 0.00198345 0.00231752 0.00257985 0.00279945 0.00298981 0.00315875 0.00331215 0.00138305 0.00187582 0.00221047 0.00247328
+          0.00269329 0.00288402 0.00305329 0.00320709 0.00122309 0.00171566 0.00205074 0.00231407 0.00253459 0.0027258 0.0028955 0.00304987
+          0.000997059 0.00148117 0.00181499 0.00207824 0.002299 0.00249055 0.00266065 0.00281574 0.000703352 0.00115165 0.00147642 0.00173614
+          0.00195523 0.00214595 0.00231566 0.00247112 0.000381658 0.000733007 0.00102236 0.00126521 0.00147468 0.00165928 0.00182484 0.00197815
+          0.000129072 0.000314535 0.000510343 0.000698089 0.00087232 0.00103266 0.00118064 0.00132152 2.02165e-05 6.63045e-05 0.000134 0.000216315
+          0.000307058 0.000401771 0.000497781 0.000596796 1.14322e-06 5.34872e-06 1.38887e-05 2.7107e-05 4.47716e-05 6.6459e-05 9.18474e-05 0.000122171
+          2.22351e-08 1.70566e-07 5.851e-07 1.36459e-06 2.55405e-06 4.17427e-06 6.26302e-06 9.05177e-06 1.80742e-10 2.87868e-09 1.34292e-08 3.66702e-08
+          7.46296e-08 1.2741e-07 1.95361e-07 2.82474e-07 9.59158e-13 3.61574e-11 2.28167e-10 7.15643e-10 1.54421e-09 2.66574e-09 3.98675e-09 5.29487e-09
+          4.97259e-15 4.17421e-13 3.43316e-12 1.21346e-11 2.74889e-11 4.76652e-11 6.88715e-11 8.35702e-11 1.38128e-16 4.97259e-15 5.04166e-14 1.95589e-13
+          4.60517e-13 8.01692e-13 1.12864e-12 1.30199e-12 0 1.38128e-16 6.90638e-16 3.31506e-15 7.87327e-15 1.35365e-14 1.86472e-14 2.11335e-14
+          0 0 1.38128e-16 0 1.38128e-16 2.76255e-16 4.14383e-16 4.14383e-16
         </DataArray>
         <DataArray type="Float32" Name="x^H2O_gas" NumberOfComponents="1" format="ascii">
-          0.00247413 0.00328858 0.00384201 0.0042765 0.00464014 0.00495534 0.00523504 0.00548917 0.00229375 0.00311034 0.00366477 0.00410009
-          0.00446446 0.0047803 0.00506057 0.00531539 0.00202868 0.00284504 0.00340027 0.00383653 0.0042018 0.00451846 0.00479949 0.00505531
-          0.001654 0.00245654 0.0030098 0.00344601 0.00381175 0.00412907 0.00441082 0.00466789 0.00116698 0.00191041 0.00244884 0.00287933
-          0.00324243 0.00355846 0.00383966 0.00409744 0.000633359 0.00121624 0.00169616 0.00209888 0.00244619 0.00275221 0.00302667 0.00328104
-          0.000214223 0.000522016 0.000846936 0.00115844 0.00144748 0.00171346 0.00195893 0.00219286 3.3555e-05 0.000110056 0.000222423 0.000359055
-          0.00050967 0.000666873 0.000826244 0.000990841 1.89753e-06 8.87781e-06 2.30538e-05 4.49969e-05 7.4323e-05 0.000110331 0.000152497 0.000202985
-          3.68839e-08 2.83036e-07 9.71057e-07 2.26495e-06 4.23963e-06 6.93001e-06 1.04011e-05 1.50544e-05 2.99533e-10 4.77505e-09 2.22822e-08 6.08543e-08
-          1.23867e-07 2.11518e-07 3.24526e-07 4.70196e-07 1.58829e-12 5.99556e-11 3.78489e-10 1.18738e-09 2.56264e-09 4.42523e-09 6.623e-09 8.81106e-09
-          8.21565e-15 6.91891e-13 5.69389e-12 2.01301e-11 4.5612e-11 7.91165e-11 1.14384e-10 1.38914e-10 0 8.21565e-15 8.34888e-14 3.24407e-13
-          7.64055e-13 1.33049e-12 1.87361e-12 2.16205e-12 0 0 1.33227e-15 5.32907e-15 1.28786e-14 2.22045e-14 3.08642e-14 3.4861e-14
-          0 0 -1.11022e-16 2.22045e-16 2.22045e-16 4.44089e-16 6.66134e-16 6.66134e-16
+          0.00239612 0.00318462 0.00372026 0.00414071 0.00449256 0.00479751 0.00506807 0.00531369 0.00222143 0.00301201 0.00354864 0.00396991
+          0.00432248 0.00462805 0.00489917 0.00514547 0.00196471 0.0027551 0.00329253 0.00371473 0.00406818 0.00437456 0.00464643 0.00489371
+          0.00160184 0.00237889 0.00291444 0.00333662 0.00369056 0.0039976 0.00427019 0.00451867 0.00113018 0.00185003 0.00237127 0.00278796
+          0.00313937 0.0034452 0.00371727 0.00396645 0.000613386 0.00117781 0.00164246 0.00203231 0.00236848 0.00266466 0.00293024 0.00317613
+          0.000207471 0.000505528 0.00082014 0.00112173 0.00140154 0.001659 0.00189656 0.00212268 3.24982e-05 0.000106582 0.000215392 0.000347688
+          0.000493513 0.000645702 0.000799958 0.000959022 1.83777e-06 8.59822e-06 2.23264e-05 4.35746e-05 7.19699e-05 0.000106831 0.00014764 0.00019638
+          3.57436e-08 2.74191e-07 9.40568e-07 2.19362e-06 4.10572e-06 6.71027e-06 1.0068e-05 1.45509e-05 2.90549e-10 4.62757e-09 2.15879e-08 5.89486e-08
+          1.1997e-07 2.04815e-07 3.14049e-07 4.54087e-07 1.54188e-12 5.81242e-11 3.66787e-10 1.15042e-09 2.48238e-09 4.28527e-09 6.40884e-09 8.51168e-09
+          7.99361e-15 6.71019e-13 5.51892e-12 1.95068e-11 4.41893e-11 7.66234e-11 1.10713e-10 1.34342e-10 2.22045e-16 7.99361e-15 8.10463e-14 3.14415e-13
+          7.40297e-13 1.28875e-12 1.81433e-12 2.09299e-12 0 2.22045e-16 1.11022e-15 5.32907e-15 1.26565e-14 2.17604e-14 2.9976e-14 3.39728e-14
+          0 0 2.22045e-16 0 2.22045e-16 4.44089e-16 6.66134e-16 6.66134e-16
         </DataArray>
         <DataArray type="Float32" Name="D^H2O_gas" NumberOfComponents="1" format="ascii">
           2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05 2.49368e-05
@@ -81,65 +107,39 @@
           2.4937e-05 2.4937e-05 2.4937e-05 2.4937e-05 2.49372e-05 2.49372e-05 2.49372e-05 2.49372e-05 2.49372e-05 2.49372e-05 2.49372e-05 2.49372e-05
           2.49374e-05 2.49374e-05 2.49374e-05 2.49374e-05 2.49374e-05 2.49374e-05 2.49374e-05 2.49374e-05
         </DataArray>
-        <DataArray type="Float32" Name="X^Air_gas" NumberOfComponents="1" format="ascii">
-          0.998459 0.997952 0.997607 0.997335 0.997108 0.996912 0.996737 0.996578 0.998572 0.998063 0.997717 0.997445
-          0.997218 0.997021 0.996846 0.996687 0.998737 0.998228 0.997882 0.99761 0.997382 0.997184 0.997009 0.996849
-          0.99897 0.99847 0.998126 0.997854 0.997625 0.997427 0.997252 0.997091 0.999274 0.998811 0.998475 0.998207
-          0.99798 0.997783 0.997608 0.997447 0.999606 0.999243 0.998944 0.998693 0.998477 0.998286 0.998115 0.997956
-          0.999867 0.999675 0.999473 0.999279 0.999099 0.998933 0.99878 0.998635 0.999979 0.999932 0.999862 0.999777
-          0.999683 0.999585 0.999486 0.999383 0.999999 0.999994 0.999986 0.999972 0.999954 0.999931 0.999905 0.999874
-          1 1 0.999999 0.999999 0.999997 0.999996 0.999994 0.999991 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1
-        </DataArray>
-        <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
-          0.997526 0.996711 0.996158 0.995723 0.99536 0.995045 0.994765 0.994511 0.997706 0.99689 0.996335 0.9959
-          0.995536 0.99522 0.994939 0.994685 0.997971 0.997155 0.9966 0.996163 0.995798 0.995482 0.995201 0.994945
-          0.998346 0.997543 0.99699 0.996554 0.996188 0.995871 0.995589 0.995332 0.998833 0.99809 0.997551 0.997121
-          0.996758 0.996442 0.99616 0.995903 0.999367 0.998784 0.998304 0.997901 0.997554 0.997248 0.996973 0.996719
-          0.999786 0.999478 0.999153 0.998842 0.998553 0.998287 0.998041 0.997807 0.999966 0.99989 0.999778 0.999641
-          0.99949 0.999333 0.999174 0.999009 0.999998 0.999991 0.999977 0.999955 0.999926 0.99989 0.999848 0.999797
-          1 1 0.999999 0.999998 0.999996 0.999993 0.99999 0.999985 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1
-        </DataArray>
-        <DataArray type="Float32" Name="velocity_liquid (m/s)" NumberOfComponents="3" format="ascii">
-          0.00578835 1.57469e-05 0 0.00620113 1.88642e-05 0 0.00610921 1.89086e-05 0 0.00601269 1.88412e-05 0
-          0.0059314 1.8749e-05 0 0.005869 1.8672e-05 0 0.00582082 1.86424e-05 0 0.00573041 1.90144e-05 0
-          0.0191001 8.89284e-06 0 0.0194189 2.02843e-05 0 0.0191317 2.06767e-05 0 0.01883 2.04216e-05 0
-          0.0185758 1.99901e-05 0 0.0183806 1.9597e-05 0 0.01823 1.93891e-05 0 0.0179557 2.13284e-05 0
-          0.0390169 -1.13593e-06 0 0.0391925 2.49863e-05 0 0.0386146 2.65327e-05 0 0.0380072 2.56598e-05 0
-          0.0374951 2.41055e-05 0 0.0371018 2.26645e-05 0 0.0367985 2.18617e-05 0 0.036268 2.86442e-05 0
-          0.068778 -1.54286e-05 0 0.0687342 3.77686e-05 0 0.067725 4.24557e-05 0 0.0666637 3.99164e-05 0
-          0.065768 3.53134e-05 0 0.0650798 3.102e-05 0 0.0645497 2.858e-05 0 0.0636756 4.7197e-05 0
-          0.113165 -3.46165e-05 0 0.112782 6.96589e-05 0 0.111136 8.2177e-05 0 0.109403 7.55278e-05 0
-          0.10794 6.33355e-05 0 0.106814 5.19145e-05 0 0.105949 4.53099e-05 0 0.10464 8.88853e-05 0
-          0.179176 -5.6722e-05 0 0.178268 0.000145624 0 0.175684 0.000176681 0 0.172966 0.0001604 0
-          0.170666 0.000130207 0 0.168895 0.000101787 0 0.167539 8.4954e-05 0 0.165724 0.000173743 0
-          0.276929 -7.07968e-05 0 0.275204 0.000321602 0 0.271252 0.000394871 0 0.267097 0.000356752 0
-          0.263574 0.000285186 0 0.260859 0.000217384 0 0.258793 0.000175718 0 0.256455 0.00032949 0
-          0.420756 -4.04749e-05 0 0.417777 0.000721519 0 0.411845 0.000887676 0 0.40562 0.000801329 0
-          0.400327 0.000636911 0 0.396244 0.00047972 0 0.393167 0.000377603 0 0.390346 0.000583985 0
-          0.630251 0.000135512 0 0.625397 0.00161463 0 0.61666 0.00197845 0 0.607517 0.0017885 0
-          0.599717 0.00142052 0 0.593688 0.001064 0 0.58921 0.000813945 0 0.585917 0.000956265 0
-          0.930517 0.000715924 0 0.922937 0.00357038 0 0.910363 0.00433653 0 0.897266 0.00393037 0
-          0.886049 0.00312804 0 0.877362 0.0023358 0 0.871031 0.00172651 0 0.867128 0.00147346 0
-          1.34957 0.00232142 0 1.33816 0.00775552 0 1.32058 0.00928304 0 1.30245 0.00843267 0
-          1.28688 0.00673193 0 1.27481 0.00501145 0 1.26621 0.00355966 0 1.26144 0.00223346 0
-          1.90783 0.00632268 0 1.89136 0.0164513 0 1.86772 0.0192443 0 1.84381 0.0174649 0
-          1.82335 0.0139689 0 1.80756 0.0103475 0 1.79655 0.0070552 0 1.79073 0.00348581 0
-          2.58722 0.0153967 0 2.56494 0.0336931 0 2.53511 0.0380765 0 2.50606 0.0342434 0
-          2.48164 0.0273051 0 2.46304 0.0200574 0 2.45036 0.0131928 0 2.44377 0.00562848 0
-          3.25094 0.0338688 0 3.22414 0.0646733 0 3.19151 0.0694232 0 3.16202 0.0609004 0
-          3.13851 0.0478768 0 3.12126 0.0346779 0 3.10981 0.0221963 0 3.10392 0.00882558 0
-          3.44724 0.064968 0 3.42491 0.10516 0 3.40468 0.103683 0 3.39082 0.0863144 0
-          3.38227 0.0655921 0 3.37724 0.0464017 0 3.37437 0.0291611 0 3.37292 0.0114325 0
-          1.96251 0.0419178 0 2.00181 0.0630285 0 2.04765 0.0594032 0 2.08784 0.0479287 0
-          2.11913 0.0356356 0 2.14177 0.0248426 0 2.15687 0.0154912 0 2.16496 0.00609638 0
+        <DataArray type="Float32" Name="velocity_gas (m/s)" NumberOfComponents="3" format="ascii">
+          0.0057881 1.51572e-05 0 0.00620058 1.82713e-05 0 0.00610858 1.83155e-05 0 0.00601203 1.82485e-05 0
+          0.00593071 1.81569e-05 0 0.00586831 1.80803e-05 0 0.00582025 1.80504e-05 0 0.0057318 1.84113e-05 0
+          0.0190993 8.31161e-06 0 0.0194172 1.96931e-05 0 0.0191298 2.00845e-05 0 0.018828 1.98295e-05 0
+          0.0185737 1.93982e-05 0 0.0183785 1.90049e-05 0 0.0182282 1.87936e-05 0 0.0179599 2.06629e-05 0
+          0.0390154 -1.69063e-06 0 0.0391893 2.44009e-05 0 0.0386109 2.59438e-05 0 0.0380033 2.50696e-05 0
+          0.037491 2.35144e-05 0 0.0370977 2.20712e-05 0 0.0367951 2.12545e-05 0 0.0362758 2.7783e-05 0
+          0.0687757 -1.59167e-05 0 0.0687291 3.71992e-05 0 0.0677189 4.18759e-05 0 0.0666571 3.93318e-05 0
+          0.0657612 3.47248e-05 0 0.065073 3.04236e-05 0 0.0645441 2.79414e-05 0 0.0636879 4.58456e-05 0
+          0.113161 -3.49579e-05 0 0.112775 6.91304e-05 0 0.111127 8.16208e-05 0 0.109394 7.49578e-05 0
+          0.107929 6.27536e-05 0 0.106804 5.13106e-05 0 0.105941 4.45946e-05 0 0.104657 8.64559e-05 0
+          0.179172 -5.67756e-05 0 0.178258 0.000145193 0 0.175673 0.000176183 0 0.172953 0.000159866 0
+          0.170652 0.000129643 0 0.168881 0.000101166 0 0.167528 8.40623e-05 0 0.165745 0.000169199 0
+          0.276925 -7.03533e-05 0 0.275193 0.000321383 0 0.271238 0.000394504 0 0.267081 0.000356304 0
+          0.263556 0.000284664 0 0.260841 0.000216728 0 0.258779 0.000174445 0 0.256477 0.000321317 0
+          0.420752 -3.92879e-05 0 0.417766 0.000721703 0 0.41183 0.000887576 0 0.405602 0.000801059 0
+          0.400308 0.000636477 0 0.396225 0.000478982 0 0.393153 0.000375557 0 0.390365 0.000570627 0
+          0.630248 0.000137658 0 0.625388 0.00161549 0 0.616647 0.00197881 0 0.6075 0.00178854 0
+          0.599698 0.00142022 0 0.59367 0.00106306 0 0.589198 0.000810476 0 0.585929 0.000937228 0
+          0.930515 0.000719129 0 0.92293 0.00357225 0 0.910352 0.00433759 0 0.897252 0.00393085 0
+          0.886033 0.00312786 0 0.877347 0.00233433 0 0.871024 0.00172077 0 0.867132 0.00145048 0
+          1.34957 0.00232556 0 1.33816 0.00775872 0 1.32057 0.00928508 0 1.30244 0.00843369 0
+          1.28686 0.00673171 0 1.2748 0.00500892 0 1.26621 0.00355108 0 1.26144 0.0022105 0
+          1.90783 0.00632734 0 1.89136 0.0164559 0 1.86772 0.0192474 0 1.8438 0.0174665 0
+          1.82334 0.0139684 0 1.80755 0.0103434 0 1.79655 0.00704441 0 1.79073 0.00346722 0
+          2.58722 0.0154011 0 2.56494 0.0336986 0 2.53511 0.0380806 0 2.50605 0.0342452 0
+          2.48164 0.027304 0 2.46304 0.020052 0 2.45036 0.013182 0 2.44377 0.00561646 0
+          3.25094 0.0338722 0 3.22414 0.0646783 0 3.19151 0.0694272 0 3.16202 0.0609021 0
+          3.13851 0.0478753 0 3.12126 0.0346728 0 3.10981 0.0221883 0 3.10392 0.00881942 0
+          3.44724 0.0649699 0 3.42491 0.105163 0 3.40468 0.103686 0 3.39083 0.0863155 0
+          3.38228 0.065591 0 3.37724 0.0463985 0 3.37438 0.029157 0 3.37292 0.0114301 0
+          1.96251 0.0419183 0 2.00181 0.0630295 0 2.04765 0.059404 0 2.08784 0.047929 0
+          2.11913 0.0356352 0 2.14177 0.0248417 0 2.15687 0.0154901 0 2.16496 0.00609581 0
         </DataArray>
         <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0