diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
index c8d3a2e6ff3adc8f3fb9eca9b45279c85b0788e9..51b065be72703049fd388daca1498453b6049739 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
@@ -709,14 +709,14 @@ class StokesDarcyCouplingDataImplementation<MDTraits, CouplingManager, enableEne
 
     using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
 
-    static constexpr bool isFicksLaw = IsFicksLaw<typename GET_PROP_TYPE(SubDomainTypeTag<stokesIdx>, MolecularDiffusionType)>();
-    static_assert(isFicksLaw == IsFicksLaw<typename GET_PROP_TYPE(SubDomainTypeTag<darcyIdx>, MolecularDiffusionType)>(),
+    static constexpr bool isFicksLaw = IsFicksLaw<GetPropType<SubDomainTypeTag<stokesIdx>, Properties::MolecularDiffusionType>>();
+    static_assert(isFicksLaw == IsFicksLaw<GetPropType<SubDomainTypeTag<darcyIdx>, Properties::MolecularDiffusionType>>(),
                   "Both submodels must use the same diffusion law.");
 
     using ReducedComponentVector = Dune::FieldVector<Scalar, numComponents-1>;
     using ReducedComponentMatrix = Dune::FieldMatrix<Scalar, numComponents-1, numComponents-1>;
 
-    using MolecularDiffusionType = typename GET_PROP_TYPE(SubDomainTypeTag<stokesIdx>, MolecularDiffusionType);
+    using MolecularDiffusionType = GetPropType<SubDomainTypeTag<stokesIdx>, Properties::MolecularDiffusionType>;
 
 public:
     using ParentType::ParentType;
@@ -918,7 +918,7 @@ protected:
      */
     Scalar diffusionCoefficientMS_(const VolumeVariables<darcyIdx>& volVars, int phaseIdx, int compIIdx, int compJIdx) const
     {
-        using EffDiffModel = typename GET_PROP_TYPE(SubDomainTypeTag<darcyIdx>, EffectiveDiffusivityModel);
+        using EffDiffModel = GetPropType<SubDomainTypeTag<darcyIdx>, Properties::EffectiveDiffusivityModel>;
         auto fluidState = volVars.fluidState();
         typename FluidSystem<darcyIdx>::ParameterCache paramCache;
         paramCache.updateAll(fluidState);
@@ -974,9 +974,9 @@ protected:
             assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx));
 
             //calculate x_inside
-            const auto xInside = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIdx);
+            const Scalar xInside = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIdx);
             //calculate outside molefraction with the respective transmissibility
-            const auto xOutside = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIdx);
+            const Scalar xOutside = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIdx);
             moleFracInside[domainICompIdx] = xInside;
             moleFracOutside[domainICompIdx] = xOutside;
         }
@@ -989,10 +989,10 @@ protected:
         for (int compIIdx = 0; compIIdx < numComponents-1; compIIdx++)
         {
             const int domainICompIIdx = couplingCompIdx(domainI, compIIdx);
-            const auto xi = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIIdx);
-            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));
+            const Scalar xi = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIIdx);
+            const Scalar avgMolarMass = volVarsI.averageMolarMass(couplingPhaseIdx(domainI));
+            const Scalar Mn = FluidSystem<i>::molarMass(numComponents-1);
+            const 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*avgMolarMass/(tin*Mn);
@@ -1005,12 +1005,12 @@ protected:
                 if (domainICompIIdx == domainICompJIdx)
                     continue;
 
-                const auto xj = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompJIdx);
-                const auto Mi = FluidSystem<i>::molarMass(domainICompIIdx);
-                const auto Mj = FluidSystem<i>::molarMass(domainICompJIdx);
-                Scalar tij = diffusionCoefficientMS_(volVarsI, couplingPhaseIdx(domainI), domainICompIIdx, domainICompJIdx);
+                const Scalar xj = volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompJIdx);
+                const Scalar Mi = FluidSystem<i>::molarMass(domainICompIIdx);
+                const Scalar Mj = FluidSystem<i>::molarMass(domainICompJIdx);
+                const Scalar tij = diffusionCoefficientMS_(volVarsI, couplingPhaseIdx(domainI), domainICompIIdx, domainICompJIdx);
                 reducedDiffusionMatrixInside[domainICompIIdx][domainICompIIdx] += xj*avgMolarMass/(tij*Mi);
-                reducedDiffusionMatrixInside[domainICompIIdx][domainICompJIdx] +=xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
+                reducedDiffusionMatrixInside[domainICompIIdx][domainICompJIdx] += xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
             }
         }
         //outside matrix
@@ -1019,10 +1019,10 @@ protected:
             const int domainJCompIIdx = couplingCompIdx(domainJ, compIIdx);
             const int domainICompIIdx = couplingCompIdx(domainI, compIIdx);
 
-            const auto xi =volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIIdx);
-            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));
+            const Scalar xi = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIIdx);
+            const Scalar avgMolarMass = volVarsJ.averageMolarMass(couplingPhaseIdx(domainJ));
+            const Scalar Mn = FluidSystem<i>::molarMass(numComponents-1);
+            const 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*avgMolarMass/(tin*Mn);
@@ -1036,17 +1036,17 @@ protected:
                 if (domainJCompJIdx == domainJCompIIdx)
                     continue;
 
-                const auto xj = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompJIdx);
-                const auto Mi = FluidSystem<i>::molarMass(domainJCompIIdx);
-                const auto Mj = FluidSystem<i>::molarMass(domainJCompJIdx);
-                Scalar tij = diffusionCoefficientMS_(volVarsJ, couplingPhaseIdx(domainJ), domainJCompIIdx, domainJCompJIdx);
+                const Scalar xj = volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompJIdx);
+                const Scalar Mi = FluidSystem<i>::molarMass(domainJCompIIdx);
+                const Scalar Mj = FluidSystem<i>::molarMass(domainJCompJIdx);
+                const Scalar tij = diffusionCoefficientMS_(volVarsJ, couplingPhaseIdx(domainJ), domainJCompIIdx, domainJCompJIdx);
                 reducedDiffusionMatrixOutside[domainICompIIdx][domainICompIIdx] += xj*avgMolarMass/(tij*Mi);
-                reducedDiffusionMatrixOutside[domainICompIIdx][domainICompJIdx] +=xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
+                reducedDiffusionMatrixOutside[domainICompIIdx][domainICompJIdx] += xi*(avgMolarMass/(tin*Mn) - avgMolarMass/(tij*Mj));
             }
         }
 
-        Scalar omegai = 1/insideDistance;
-        Scalar omegaj = 1/outsideDistance;
+        const Scalar omegai = 1/insideDistance;
+        const Scalar omegaj = 1/outsideDistance;
 
         reducedDiffusionMatrixInside.invert();
         reducedDiffusionMatrixInside *= omegai*volVarsI.density(couplingPhaseIdx(domainI));
@@ -1095,25 +1095,22 @@ protected:
     {
         NumEqVector diffusiveFlux(0.0);
 
-        const auto rhoInside = (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) ? volVarsI.density(couplingPhaseIdx(domainI)) :  volVarsI.molarDensity(couplingPhaseIdx(domainI));
-
-        const auto rhoOutside = (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) ? volVarsJ.density(couplingPhaseIdx(domainJ)) :  volVarsJ.molarDensity(couplingPhaseIdx(domainJ));
-
+        const Scalar rhoInside = massOrMolarDensity(volVarsI, MolecularDiffusionType::referenceSystemFormulation(), couplingPhaseIdx(domainI));
+        const Scalar rhoOutside = massOrMolarDensity(volVarsJ, MolecularDiffusionType::referenceSystemFormulation(), couplingPhaseIdx(domainJ));
         const Scalar avgDensity = 0.5 * rhoInside + 0.5 * rhoOutside;
 
         const Scalar insideDistance = this->getDistance_(scvI, scvfI);
         const Scalar outsideDistance = this->getDistance_(scvJ, scvfI);
 
-        for(int compIdx = 1; compIdx < numComponents; ++compIdx)
+        for (int compIdx = 1; compIdx < numComponents; ++compIdx)
         {
             const int domainICompIdx = couplingCompIdx(domainI, compIdx);
             const int domainJCompIdx = couplingCompIdx(domainJ, compIdx);
 
             assert(FluidSystem<i>::componentName(domainICompIdx) == FluidSystem<j>::componentName(domainJCompIdx));
 
-            const auto massOrMoleFractionInside = (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) ?volVarsI.massFraction(couplingPhaseIdx(domainI), domainICompIdx) :  volVarsI.moleFraction(couplingPhaseIdx(domainI), domainICompIdx);
-
-            const auto massOrMoleFractionOutside = (MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged) ?volVarsJ.massFraction(couplingPhaseIdx(domainJ), domainJCompIdx):  volVarsJ.moleFraction(couplingPhaseIdx(domainJ), domainJCompIdx);
+            const Scalar massOrMoleFractionInside = massOrMoleFraction(volVarsI, MolecularDiffusionType::referenceSystemFormulation(), couplingPhaseIdx(domainI), domainICompIdx);
+            const Scalar massOrMoleFractionOutside = massOrMoleFraction(volVarsJ, MolecularDiffusionType::referenceSystemFormulation(), couplingPhaseIdx(domainJ), domainJCompIdx);
 
             const Scalar deltaMassOrMoleFrac = massOrMoleFractionOutside - massOrMoleFractionInside;
             const Scalar tij = this->transmissibility_(domainI,