diff --git a/dumux/flux/box/fickslaw.hh b/dumux/flux/box/fickslaw.hh
index 3893f733b66b94f3b0bffd87181cf19bee619e87..a0abc1452853cae9d369dcd3353ad831ee42a1c6 100644
--- a/dumux/flux/box/fickslaw.hh
+++ b/dumux/flux/box/fickslaw.hh
@@ -20,7 +20,7 @@
  * \file
  * \ingroup BoxFlux
  * \brief This file contains the data which is required to calculate
- *        diffusive mass fluxes due to molecular diffusion with Fick's law.
+ *        diffusive fluxes due to molecular diffusion with Fick's law.
  */
 #ifndef DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
 #define DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
@@ -71,7 +71,8 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::box, referenceSystem
 
 public:
     //return the reference system
-    static constexpr ReferenceSystemFormulation referenceSystemFormulation() { return referenceSystem; }
+    static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+    { return referenceSystem; }
 
     static ComponentFluxVector flux(const Problem& problem,
                                     const Element& element,
diff --git a/dumux/flux/box/maxwellstefanslaw.hh b/dumux/flux/box/maxwellstefanslaw.hh
index ed33bdaa00158dc51c21dceb121531102126b1b6..ceb50bda2cc61b74dc1acd05843c5be3f626b63d 100644
--- a/dumux/flux/box/maxwellstefanslaw.hh
+++ b/dumux/flux/box/maxwellstefanslaw.hh
@@ -70,7 +70,8 @@ class MaxwellStefansLawImplementation<TypeTag, DiscretizationMethod::box, refere
 
 public:
     //return the reference system
-    static constexpr ReferenceSystemFormulation referenceSystemFormulation() { return referenceSystem; }
+    static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+    { return referenceSystem; }
 
     static ComponentFluxVector flux(const Problem& problem,
                                     const Element& element,
@@ -93,7 +94,7 @@ public:
         const auto& shapeValues = fluxVarsCache.shapeValues();
 
         Scalar rho(0.0);
-        Scalar Mavg(0.0);
+        Scalar avgMolarMass(0.0);
         for (auto&& scv : scvs(fvGeometry))
         {
             const auto& volVars = elemVolVars[scv];
@@ -101,7 +102,7 @@ public:
             // density interpolation
             rho +=  volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0];
             //average molar mass interpolation
-            Mavg += volVars.averageMolarMass(phaseIdx)*shapeValues[scv.indexInElement()][0];
+            avgMolarMass += volVars.averageMolarMass(phaseIdx)*shapeValues[scv.indexInElement()][0];
             //interpolate the mole fraction for the diffusion matrix
             for (int compIdx = 0; compIdx < numComponents; compIdx++)
             {
@@ -109,7 +110,7 @@ public:
             }
         }
 
-        reducedDiffusionMatrix = setupMSMatrix_(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, moleFrac, Mavg);
+        reducedDiffusionMatrix = setupMSMatrix_(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, moleFrac, avgMolarMass);
 
         for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
         {
@@ -144,7 +145,7 @@ private:
                                                  const SubControlVolumeFace& scvf,
                                                  const int phaseIdx,
                                                  const ComponentFluxVector moleFrac,
-                                                 const Scalar Mavg)
+                                                 const Scalar avgMolarMass)
     {
         ReducedComponentMatrix reducedDiffusionMatrix(0.0);
 
@@ -178,7 +179,7 @@ private:
             const auto tin = problem.spatialParams().harmonicMean(tinInside, tinOutside, scvf.unitOuterNormal());
 
             //begin the entrys of the diffusion matrix of the diagonal
-            reducedDiffusionMatrix[compIIdx][compIIdx] += xi*Mavg/(tin*Mn);
+            reducedDiffusionMatrix[compIIdx][compIIdx] += xi*avgMolarMass/(tin*Mn);
 
             // now set the rest of the entries (off-diagonal and additional entries for diagonal)
             for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
@@ -203,9 +204,9 @@ private:
                 // the resulting averaged diffusion tensor
                 const auto tij = problem.spatialParams().harmonicMean(tijInside, tijOutside, scvf.unitOuterNormal());
 
-                reducedDiffusionMatrix[compIIdx][compIIdx] += xj*Mavg/(tij*Mi);
+                reducedDiffusionMatrix[compIIdx][compIIdx] += xj*avgMolarMass/(tij*Mi);
                 if (compJIdx < numComponents-1)
-                    reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(Mavg/(tin*Mn) - Mavg/(tij*Mj));
+                    reducedDiffusionMatrix[compIIdx][compJIdx] +=xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
             }
         }
         return reducedDiffusionMatrix;
diff --git a/dumux/flux/ccmpfa/fickslaw.hh b/dumux/flux/ccmpfa/fickslaw.hh
index 26a34ed3d56a70cca448f4e3b6e56c08a43bf82c..2d8cc2d2204ade8ebef337dbaf976156d456c048 100644
--- a/dumux/flux/ccmpfa/fickslaw.hh
+++ b/dumux/flux/ccmpfa/fickslaw.hh
@@ -161,7 +161,8 @@ public:
     // state the discretization method this implementation belongs to
     static const DiscretizationMethod discMethod = DiscretizationMethod::ccmpfa;
     //return the reference system
-    static constexpr ReferenceSystemFormulation referenceSystemFormulation() { return referenceSystem; }
+    static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+    { return referenceSystem; }
     // state the type for the corresponding cache and its filler
     using Cache = MpfaFicksLawCache;
 
diff --git a/dumux/flux/cctpfa/fickslaw.hh b/dumux/flux/cctpfa/fickslaw.hh
index d065233e6b930af57f67b9c76d272820bb283d9b..c63c21eebf0f0fe063eb2d762e04e2ed944bb2a3 100644
--- a/dumux/flux/cctpfa/fickslaw.hh
+++ b/dumux/flux/cctpfa/fickslaw.hh
@@ -117,7 +117,8 @@ public:
     //! state the discretization method this implementation belongs to
     static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
     //! Return the reference system
-    static constexpr ReferenceSystemFormulation referenceSystemFormulation() { return referenceSystem; }
+    static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+    { return referenceSystem; }
 
     //! state the type for the corresponding cache and its filler
     using Cache = TpfaFicksLawCache;
diff --git a/dumux/flux/cctpfa/maxwellstefanslaw.hh b/dumux/flux/cctpfa/maxwellstefanslaw.hh
index fd3c85dad4f35224661621907d8b5ff8420b6d2f..c0682083e09212089ae5ff367e9c3dd9abcf2716 100644
--- a/dumux/flux/cctpfa/maxwellstefanslaw.hh
+++ b/dumux/flux/cctpfa/maxwellstefanslaw.hh
@@ -76,7 +76,8 @@ public:
     // state the discretization method this implementation belongs to
     static const DiscretizationMethod discMethod = DiscretizationMethod::cctpfa;
     //return the reference system
-    static constexpr ReferenceSystemFormulation referenceSystemFormulation() { return referenceSystem; }
+    static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+    { return referenceSystem; }
 
     //! state the type for the corresponding cache and its filler
     //! We don't cache anything for this law
@@ -217,7 +218,7 @@ private:
         if(Dune::FloatCmp::eq<Scalar>(volVars.saturation(phaseIdx), 0))
             return reducedDiffusionMatrix;
 
-        const auto Mavg = volVars.averageMolarMass(phaseIdx);
+        const auto avgMolarMass = volVars.averageMolarMass(phaseIdx);
 
         for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
         {
@@ -228,7 +229,7 @@ private:
             tin = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tin);
 
             // set the entries of the diffusion matrix of the diagonal
-            reducedDiffusionMatrix[compIIdx][compIIdx] += xi*Mavg/(tin*Mn);
+            reducedDiffusionMatrix[compIIdx][compIIdx] += xi*avgMolarMass/(tin*Mn);
 
             // now set the rest of the entries (off-diagonal and additional entries for diagonal)
             for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
@@ -242,9 +243,9 @@ private:
                 const auto Mj = FluidSystem::molarMass(compJIdx);
                 Scalar tij = getDiffusionCoefficient(phaseIdx, compIIdx, compJIdx, problem, element, volVars, scv);
                 tij = EffDiffModel::effectiveDiffusivity(volVars.porosity(), volVars.saturation(phaseIdx), tij);
-                reducedDiffusionMatrix[compIIdx][compIIdx] += xj*Mavg/(tij*Mi);
+                reducedDiffusionMatrix[compIIdx][compIIdx] += xj*avgMolarMass/(tij*Mi);
                 if (compJIdx < numComponents-1)
-                    reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(Mavg/(tin*Mn) - Mavg/(tij*Mj));
+                    reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
             }
         }
         return reducedDiffusionMatrix;
diff --git a/dumux/flux/fickslaw.hh b/dumux/flux/fickslaw.hh
index d3afa33925fbbbbe4ff8ea64df251959af0e93af..ad5aa09ed9b09e14fcdbd094fae377b7326bd4b4 100644
--- a/dumux/flux/fickslaw.hh
+++ b/dumux/flux/fickslaw.hh
@@ -28,7 +28,7 @@
 
 #include <dumux/common/properties.hh>
 #include <dumux/discretization/method.hh>
-#include <dumux/flux/referencesystem.hh>
+#include <dumux/flux/referencesystemformulation.hh>
 
 namespace Dumux {
 
diff --git a/dumux/flux/maxwellstefanslaw.hh b/dumux/flux/maxwellstefanslaw.hh
index 288f6688cede38c734ec26130274bf3b0d820818..6bf9966bf22e40391c1d3dec5d8b236183b05203 100644
--- a/dumux/flux/maxwellstefanslaw.hh
+++ b/dumux/flux/maxwellstefanslaw.hh
@@ -27,7 +27,7 @@
 
 #include <dumux/common/properties.hh>
 #include <dumux/discretization/method.hh>
-#include <dumux/flux/referencesystem.hh>
+#include <dumux/flux/referencesystemformulation.hh>
 
 namespace Dumux {
 
diff --git a/dumux/flux/referencesystem.hh b/dumux/flux/referencesystemformulation.hh
similarity index 54%
rename from dumux/flux/referencesystem.hh
rename to dumux/flux/referencesystemformulation.hh
index d37f6af8b638fa56c31d6270be8a5de854ce5ded..8ca6c0697c3d4542b18d40a99dac5ad9df5b3705 100644
--- a/dumux/flux/referencesystem.hh
+++ b/dumux/flux/referencesystemformulation.hh
@@ -22,22 +22,22 @@
  * \brief The reference frameworks and formulations available for splitting total fluxes into a advective and diffusive part
  */
 
-#ifndef DUMUX_FLUX_REFERENCESYSTEM_HH
-#define DUMUX_FLUX_REFERENCESYSTEM_HH
+#ifndef DUMUX_FLUX_REFERENCESYSTEMFORMULATION_HH
+#define DUMUX_FLUX_REFERENCESYSTEMFORMULATION_HH
 
 namespace Dumux {
 
-    /*!
-     * \brief The formulations available for Fick's law related to the reference system
-     * \ingroup Flux
-     * \note The total flux of a component can be split into an advective and a diffusive part. The advective part is in our framework determined by a momentum balance. Standard momentum balances, e.g. Navier-Stokes and Darcy's law give back mass-averaged velocites (see Multicomponent Mass Transfer, R. Taylor u. R. Krishna. J.), therefore as default we use the formulation of Fick's law which matches to these velocities (mass averaged formulation).
-     *
-     * This means that the diffusive fluxes are calculated with the mass fraction gradients and the unit of the fluxes is in kg/s. It is also possible to use a molar averaged reference system, which can be benefitial e.g. when it is known that the molar averaged advective velocity would be zero. When using a molar averaged reference velocity Fick's law is formulated with mole fraction gradients and the unit of the flux is mol/s. This means that depending on the reference system the units of the fluxes need to be adapted to be used in mass or mole balances.
-     */
-    enum class ReferenceSystemFormulation
-    {
-        massAveraged, molarAveraged
-    };
+/*!
+ * \brief The formulations available for Fick's law related to the reference system
+ * \ingroup Flux
+ * \note The total flux of a component can be split into an advective and a diffusive part. The advective part is in our framework determined by a momentum balance. Standard momentum balances, e.g. Navier-Stokes and Darcy's law give back mass-averaged velocites (see Multicomponent Mass Transfer, R. Taylor u. R. Krishna. J.), therefore as default we use the formulation of Fick's law which matches to these velocities (mass averaged formulation).
+ *
+ * This means that the diffusive fluxes are calculated with the mass fraction gradients and the unit of the fluxes is in kg/s. It is also possible to use a molar averaged reference system, which can be benefitial e.g. when it is known that the molar averaged advective velocity would be zero. When using a molar averaged reference velocity Fick's law is formulated with mole fraction gradients and the unit of the flux is mol/s. This means that depending on the reference system the units of the fluxes need to be adapted to be used in mass or mole balances.
+*/
+enum class ReferenceSystemFormulation
+{
+    massAveraged, molarAveraged
+};
 
 } // end namespace Dumux
 
diff --git a/dumux/flux/staggered/freeflow/fickslaw.hh b/dumux/flux/staggered/freeflow/fickslaw.hh
index dbd01b1071bf7c37570d298cbf87586d809a01a5..68b33982fac1e340388802df8d6976a7b6573e59 100644
--- a/dumux/flux/staggered/freeflow/fickslaw.hh
+++ b/dumux/flux/staggered/freeflow/fickslaw.hh
@@ -68,7 +68,8 @@ public:
     // state the discretization method this implementation belongs to
     static const DiscretizationMethod discMethod = DiscretizationMethod::staggered;
     //return the reference system
-    static constexpr ReferenceSystemFormulation referenceSystemFormulation() { return referenceSystem; }
+    static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+    { return referenceSystem; }
 
     //! state the type for the corresponding cache
     //! We don't cache anything for this law
diff --git a/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh b/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh
index b499a9dc5d7a303bbf9cab8e2054bc0b20ccab73..d6534fa82b3b25ad00297a2844520d4fe774ab56 100644
--- a/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh
+++ b/dumux/flux/staggered/freeflow/maxwellstefanslaw.hh
@@ -74,7 +74,8 @@ public:
     // state the discretization method this implementation belongs to
     static const DiscretizationMethod discMethod = DiscretizationMethod::staggered;
     //return the reference system
-    static constexpr ReferenceSystemFormulation referenceSystemFormulation() { return referenceSystem; }
+    static constexpr ReferenceSystemFormulation referenceSystemFormulation()
+    { return referenceSystem; }
 
     //! state the type for the corresponding cache and its filler
     //! We don't cache anything for this law
@@ -212,12 +213,12 @@ private:
         for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
         {
             const auto xi = volVars.moleFraction(compIIdx);
-            const auto Mavg = volVars.averageMolarMass(0);
+            const auto avgMolarMass = volVars.averageMolarMass(0);
             const auto Mn = FluidSystem::molarMass(numComponents-1);
             const Scalar tin = volVars.effectiveDiffusivity(compIIdx, numComponents-1);
 
             // set the entries of the diffusion matrix of the diagonal
-            reducedDiffusionMatrix[compIIdx][compIIdx] +=  xi*Mavg/(tin*Mn);
+            reducedDiffusionMatrix[compIIdx][compIIdx] +=  xi*avgMolarMass/(tin*Mn);
 
             for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
             {
@@ -229,9 +230,9 @@ private:
                 const auto Mi = FluidSystem::molarMass(compIIdx);
                 const auto Mj = FluidSystem::molarMass(compJIdx);
                 const Scalar tij = volVars.effectiveDiffusivity(compIIdx, compJIdx);
-                reducedDiffusionMatrix[compIIdx][compIIdx] +=  xj*Mavg/(tij*Mi);
+                reducedDiffusionMatrix[compIIdx][compIIdx] +=  xj*avgMolarMass/(tij*Mi);
                 if (compJIdx < numComponents-1)
-                    reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(Mavg/(tin*Mn) - Mavg/(tij*Mj));
+                    reducedDiffusionMatrix[compIIdx][compJIdx] += xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
             }
         }
         return reducedDiffusionMatrix;
diff --git a/dumux/freeflow/compositional/staggered/fluxvariables.hh b/dumux/freeflow/compositional/staggered/fluxvariables.hh
index de283bab6fdc6f9fff9e24100201bdcedd8138e4..3eab6b5ab0b9ef06d1e1c349a7b8b5bc7c64bc32 100644
--- a/dumux/freeflow/compositional/staggered/fluxvariables.hh
+++ b/dumux/freeflow/compositional/staggered/fluxvariables.hh
@@ -75,6 +75,8 @@ public:
 
         const auto diffusiveFluxes = MolecularDiffusionType::flux(problem, element, fvGeometry, elemVolVars, scvf);
 
+        auto referenceSystemFormulation = MolecularDiffusionType::referenceSystemFormulation();
+
         for (int compIdx = 0; compIdx < numComponents; ++compIdx)
         {
             auto upwindTerm = [compIdx](const auto& volVars)
@@ -87,12 +89,14 @@ public:
             flux[compIdx] = ParentType::advectiveFluxForCellCenter(problem, elemVolVars, elemFaceVars, scvf, upwindTerm);
 
             //check for the reference system and adapt units of the diffusive flux accordingly.
-            if (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+            if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
             {
                 flux[compIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
             }
-            else
+            else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
                 flux[compIdx] += useMoles ? diffusiveFluxes[compIdx] : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+            else
+                DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
         }
 
         // in case one balance is substituted by the total mass balance
diff --git a/dumux/freeflow/nonisothermal/localresidual.hh b/dumux/freeflow/nonisothermal/localresidual.hh
index 0a5475a3a758fa19e6ac90e199958dd8a9eb7513..abb3676f68294f2d0ed605f5558157065cc08f8c 100644
--- a/dumux/freeflow/nonisothermal/localresidual.hh
+++ b/dumux/freeflow/nonisothermal/localresidual.hh
@@ -25,7 +25,7 @@
 #define DUMUX_FREE_FLOW_ENERGY_LOCAL_RESIDUAL_HH
 
 #include <dumux/discretization/method.hh>
-#include <dumux/flux/referencesystem.hh>
+#include <dumux/flux/referencesystemformulation.hh>
 
 namespace Dumux {
 
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
index 1af4a2673ba76707cd4ef9af1bd170fcc6c3dbbc..c8d3a2e6ff3adc8f3fb9eca9b45279c85b0788e9 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
@@ -30,7 +30,7 @@
 #include <dumux/common/properties.hh>
 #include <dumux/common/math.hh>
 #include <dumux/discretization/method.hh>
-#include <dumux/flux/referencesystem.hh>
+#include <dumux/flux/referencesystemformulation.hh>
 #include <dumux/multidomain/couplingmanager.hh>
 #include <dumux/common/deprecated.hh>
 
@@ -990,12 +990,12 @@ protected:
         {
             const int domainICompIIdx = couplingCompIdx(domainI, compIIdx);
             const auto xi = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIIdx);
-            const auto Mavg = volVarsI.fluidState().averageMolarMass(couplingPhaseIdx(domainI));
+            const auto avgMolarMass = volVarsI.fluidState().averageMolarMass(couplingPhaseIdx(domainI));
             const auto Mn = FluidSystem<i>::molarMass(numComponents-1);
             Scalar tin = diffusionCoefficientMS_(volVarsI, couplingPhaseIdx(domainI), domainICompIIdx, couplingCompIdx(domainI, numComponents-1));
 
             // set the entries of the diffusion matrix of the diagonal
-            reducedDiffusionMatrixInside[domainICompIIdx][domainICompIIdx] += xi*Mavg/(tin*Mn);
+            reducedDiffusionMatrixInside[domainICompIIdx][domainICompIIdx] += xi*avgMolarMass/(tin*Mn);
 
             for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
             {
@@ -1009,8 +1009,8 @@ protected:
                 const auto Mi = FluidSystem<i>::molarMass(domainICompIIdx);
                 const auto Mj = FluidSystem<i>::molarMass(domainICompJIdx);
                 Scalar tij = diffusionCoefficientMS_(volVarsI, couplingPhaseIdx(domainI), domainICompIIdx, domainICompJIdx);
-                reducedDiffusionMatrixInside[domainICompIIdx][domainICompIIdx] += xj*Mavg/(tij*Mi);
-                reducedDiffusionMatrixInside[domainICompIIdx][domainICompJIdx] +=xi*(Mavg/(tin*Mn) - Mavg/(tij*Mj));
+                reducedDiffusionMatrixInside[domainICompIIdx][domainICompIIdx] += xj*avgMolarMass/(tij*Mi);
+                reducedDiffusionMatrixInside[domainICompIIdx][domainICompJIdx] +=xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
             }
         }
         //outside matrix
@@ -1020,12 +1020,12 @@ protected:
             const int domainICompIIdx = couplingCompIdx(domainI, compIIdx);
 
             const auto xi =volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIIdx);
-            const auto Mavg = volVarsJ.fluidState().averageMolarMass(couplingPhaseIdx(domainJ));
+            const auto avgMolarMass = volVarsJ.fluidState().averageMolarMass(couplingPhaseIdx(domainJ));
             const auto Mn = FluidSystem<i>::molarMass(numComponents-1);
             Scalar tin =diffusionCoefficientMS_(volVarsJ, couplingPhaseIdx(domainJ), domainJCompIIdx, couplingCompIdx(domainJ, numComponents-1));
 
             // set the entries of the diffusion matrix of the diagonal
-            reducedDiffusionMatrixOutside[domainICompIIdx][domainICompIIdx] +=  xi*Mavg/(tin*Mn);
+            reducedDiffusionMatrixOutside[domainICompIIdx][domainICompIIdx] +=  xi*avgMolarMass/(tin*Mn);
 
             for (int compJIdx = 0; compJIdx < numComponents; compJIdx++)
             {
@@ -1040,8 +1040,8 @@ protected:
                 const auto Mi = FluidSystem<i>::molarMass(domainJCompIIdx);
                 const auto Mj = FluidSystem<i>::molarMass(domainJCompJIdx);
                 Scalar tij = diffusionCoefficientMS_(volVarsJ, couplingPhaseIdx(domainJ), domainJCompIIdx, domainJCompJIdx);
-                reducedDiffusionMatrixOutside[domainICompIIdx][domainICompIIdx] += xj*Mavg/(tij*Mi);
-                reducedDiffusionMatrixOutside[domainICompIIdx][domainICompJIdx] +=xi*(Mavg/(tin*Mn) - Mavg/(tij*Mj));
+                reducedDiffusionMatrixOutside[domainICompIIdx][domainICompIIdx] += xj*avgMolarMass/(tij*Mi);
+                reducedDiffusionMatrixOutside[domainICompIIdx][domainICompJIdx] +=xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
             }
         }
 
diff --git a/dumux/porousmediumflow/3p3c/localresidual.hh b/dumux/porousmediumflow/3p3c/localresidual.hh
index a89071675144c33e9fbed6b8825b6f13e165526d..750d25e0731d3a300b975e2cae2776c2e2293847 100644
--- a/dumux/porousmediumflow/3p3c/localresidual.hh
+++ b/dumux/porousmediumflow/3p3c/localresidual.hh
@@ -175,11 +175,13 @@ public:
             jGW += diffusionFluxesWPhase[gCompIdx]/FluidSystem::molarMass(gCompIdx);
             jNW += diffusionFluxesWPhase[nCompIdx]/FluidSystem::molarMass(nCompIdx);
         }
-        else
+        else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
         {
             jGW += diffusionFluxesWPhase[gCompIdx];
             jNW += diffusionFluxesWPhase[nCompIdx];
         }
+        else
+            DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
 
         jWW += -(jGW+jNW);
 
diff --git a/dumux/porousmediumflow/3pwateroil/localresidual.hh b/dumux/porousmediumflow/3pwateroil/localresidual.hh
index 42cd8efbf1bd4f2dbd718d6af99fa35eb4baa664..c3d360ffc9271f8ea6adc70a9682dbe82a58fd62 100644
--- a/dumux/porousmediumflow/3pwateroil/localresidual.hh
+++ b/dumux/porousmediumflow/3pwateroil/localresidual.hh
@@ -174,14 +174,18 @@ public:
 
         const auto diffusionFluxesWPhase = fluxVars.molecularDiffusionFlux(wPhaseIdx);
 
+        auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
+
         Scalar jNW = 0.0;
         Scalar jWW = 0.0;
         // diffusive fluxes
         //check for the reference system and adapt units of the diffusive flux accordingly.
-        if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+        if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
             jNW +=  diffusionFluxesWPhase[nCompIdx]/FluidSystem::molarMass(nCompIdx);
-        else
+        else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
             jNW +=  diffusionFluxesWPhase[nCompIdx];
+        else
+            DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
 
         jWW += -(jNW);
 
@@ -190,7 +194,7 @@ public:
         Scalar jWG = 0.0;
         Scalar jNG = 0.0;
 
-        if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+        if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
             jWG += diffusionFluxesGPhase[wCompIdx]/FluidSystem::molarMass(wCompIdx);
         else
             jWG += diffusionFluxesGPhase[wCompIdx];
@@ -202,7 +206,7 @@ public:
         Scalar jWN = 0.0;
         Scalar jNN = 0.0;
 
-        if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+        if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
             jWN += diffusionFluxesNPhase[wCompIdx]/FluidSystem::molarMass(wCompIdx);
         else
             jWN += diffusionFluxesNPhase[wCompIdx];
diff --git a/dumux/porousmediumflow/compositional/localresidual.hh b/dumux/porousmediumflow/compositional/localresidual.hh
index f878772ead911e7e8fcdb7d40a870d038b7b0297..9ad1fa80fc7cbdcdb581e09e30a1fc613e78169e 100644
--- a/dumux/porousmediumflow/compositional/localresidual.hh
+++ b/dumux/porousmediumflow/compositional/localresidual.hh
@@ -176,9 +176,11 @@ public:
                     if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
                         flux[eqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx)
                                             : diffusiveFluxes[compIdx];
-                    else
+                    else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
                         flux[eqIdx] += useMoles ? diffusiveFluxes[compIdx]
                                             : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+                    else
+                        DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
                 }
             }
 
@@ -196,9 +198,11 @@ public:
                     //check for the reference system and adapt units of the diffusive flux accordingly.
                     if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
                         flux[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx) : diffusiveFluxes[compIdx];
-                    else
+                    else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
                         flux[replaceCompEqIdx] += useMoles ? diffusiveFluxes[compIdx]
                                                 : diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+                    else
+                        DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
                 }
             }
 
diff --git a/dumux/porousmediumflow/fluxvariablescachefiller.hh b/dumux/porousmediumflow/fluxvariablescachefiller.hh
index 22f6896209e3d242468cdfd764193024d3bacac4..40ec6db696d9990615d0d089e76994f38c79391e 100644
--- a/dumux/porousmediumflow/fluxvariablescachefiller.hh
+++ b/dumux/porousmediumflow/fluxvariablescachefiller.hh
@@ -28,7 +28,7 @@
 #include <dumux/common/parameters.hh>
 
 #include <dumux/discretization/method.hh>
-#include <dumux/flux/referencesystem.hh>
+#include <dumux/flux/referencesystemformulation.hh>
 #include <dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh>
 
 namespace Dumux {
diff --git a/dumux/porousmediumflow/tracer/localresidual.hh b/dumux/porousmediumflow/tracer/localresidual.hh
index 538eeaf18dfe759cdc2bd7ee085e4e43287321df..a56ec4e0512b011a087f141d2bdc4f64800ae13f 100644
--- a/dumux/porousmediumflow/tracer/localresidual.hh
+++ b/dumux/porousmediumflow/tracer/localresidual.hh
@@ -126,6 +126,7 @@ public:
         NumEqVector flux(0.0);
         const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
 
+        auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
         // formulation with mole balances
         if (useMoles)
         {
@@ -138,10 +139,12 @@ public:
                 // advective fluxes
                 flux[compIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
                 // diffusive fluxes
-                if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+                if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
                     flux[compIdx] += diffusiveFluxes[compIdx]/FluidSystem::molarMass(compIdx);
-                else
+                else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
                     flux[compIdx] += diffusiveFluxes[compIdx];
+                else
+                    DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
             }
         }
         // formulation with mass balances
@@ -156,10 +159,12 @@ public:
                 // advective fluxes
                 flux[compIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
                 // diffusive fluxes
-                if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+                if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
                     flux[compIdx] += diffusiveFluxes[compIdx];
-                else
+                else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
                     flux[compIdx] += diffusiveFluxes[compIdx]*FluidSystem::molarMass(compIdx);
+                else
+                    DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
             }
         }
 
@@ -260,11 +265,13 @@ public:
                 diffDeriv = useMoles ? massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx)/FluidSystem::molarMass(compIdx)
                                             : massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx);
             }
-            else
+            else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
             {
                 diffDeriv = useMoles ? massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx)
                                             : massOrMolarDensity*fluxCache.diffusionTij(phaseIdx,         compIdx)*FluidSystem::molarMass(compIdx);
             }
+            else
+                DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
 
             derivativeMatrices[scvf.insideScvIdx()][compIdx][compIdx] += (advDerivII + diffDeriv);
             derivativeMatrices[scvf.outsideScvIdx()][compIdx][compIdx] += (advDerivIJ - diffDeriv);
@@ -323,9 +330,11 @@ public:
                 if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
                     diffDeriv += useMoles ? ti[compIdx][scv.indexInElement()]/FluidSystem::molarMass(compIdx)
                                         : ti[compIdx][scv.indexInElement()];
-                else
+                else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
                     diffDeriv += useMoles ? ti[compIdx][scv.indexInElement()]
                                             : ti[compIdx][scv.indexInElement()]*FluidSystem::molarMass(compIdx);
+                else
+                    DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
                 A[insideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] += diffDeriv;
                 A[outsideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] -= diffDeriv;
             }
diff --git a/test/porousmediumflow/tracer/multicomp/problem.hh b/test/porousmediumflow/tracer/multicomp/problem.hh
index af9e1e49ddd18ee67b8dd462b13c06b0bb3378d6..90acc969b74db95d60c9bdc217b9b413ca782102 100644
--- a/test/porousmediumflow/tracer/multicomp/problem.hh
+++ b/test/porousmediumflow/tracer/multicomp/problem.hh
@@ -97,18 +97,18 @@ public:
     //! The number of components
     static constexpr int numComponents = 3;
 
-    static constexpr int CompOneIdx = 0;//first major component
-    static constexpr int CompTwoIdx = 1;//second major component
-    static constexpr int CompThreeIdx = 2;//secondary component
+    static constexpr int compOneIdx = 0;//first major component
+    static constexpr int compTwoIdx = 1;//second major component
+    static constexpr int compThreeIdx = 2;//secondary component
 
     //! Human readable component name (index compIdx) (for vtk output)
     static std::string componentName(int compIdx)
     {
         switch (compIdx)
         {
-        case CompOneIdx: return "CompOne";
-        case CompTwoIdx: return "CompTwo";
-        case CompThreeIdx:return "CompThree";
+        case compOneIdx: return "CompOne";
+        case compTwoIdx: return "CompTwo";
+        case compThreeIdx:return "CompThree";
         }
         DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
     }
@@ -128,11 +128,11 @@ public:
                                              const Element& element,
                                              const SubControlVolume& scv)
     {
-      if (compIdx == CompOneIdx)
+      if (compIdx == compOneIdx)
           return 0;
-      if (compIdx == CompTwoIdx)
+      if (compIdx == compTwoIdx)
           return 83.3e-6;
-      if (compIdx == CompThreeIdx)
+      if (compIdx == compThreeIdx)
           return 68.0e-6;
        DUNE_THROW(Dune::InvalidStateException,
                        "Binary diffusion coefficient of component "
@@ -153,11 +153,11 @@ public:
             swap(compIIdx, compJIdx);
         }
 
-        if (compIIdx == CompOneIdx && compJIdx == CompTwoIdx)
+        if (compIIdx == compOneIdx && compJIdx == compTwoIdx)
             return 83.3e-6;
-        if (compIIdx == CompOneIdx && compJIdx == CompThreeIdx)
+        if (compIIdx == compOneIdx && compJIdx == compThreeIdx)
             return 68.0e-6;
-        if (compIIdx == CompTwoIdx && compJIdx == CompThreeIdx)
+        if (compIIdx == compTwoIdx && compJIdx == compThreeIdx)
             return 16.8e-6;
         DUNE_THROW(Dune::InvalidStateException,
                        "Binary diffusion coefficient of components "
@@ -385,15 +385,15 @@ public:
         PrimaryVariables initialValues(0.0);
         if (globalPos[0] < 0.5)
         {
-           initialValues[FluidSystem::CompOneIdx] = 0.0;
-           initialValues[FluidSystem::CompTwoIdx] = 0.50086;
-           initialValues[FluidSystem::CompThreeIdx] = 0.49914;
+           initialValues[FluidSystem::compOneIdx] = 0.0;
+           initialValues[FluidSystem::compTwoIdx] = 0.50086;
+           initialValues[FluidSystem::compThreeIdx] = 0.49914;
         }
         else
         {
-           initialValues[FluidSystem::CompOneIdx] = 0.50121;
-           initialValues[FluidSystem::CompTwoIdx] = 0.49879;
-           initialValues[FluidSystem::CompThreeIdx] = 0.0;
+           initialValues[FluidSystem::compOneIdx] = 0.50121;
+           initialValues[FluidSystem::compTwoIdx] = 0.49879;
+           initialValues[FluidSystem::compThreeIdx] = 0.0;
         }
         return initialValues;
     }