diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh
index a072788125373e4451b081c9a93b79e0ba8a8b88..56b4ee040a5e092d595dcc7a9e41d54bd030748f 100644
--- a/dumux/freeflow/stokes/stokesfluxvariables.hh
+++ b/dumux/freeflow/stokes/stokesfluxvariables.hh
@@ -297,6 +297,12 @@ public:
     const DimMatrix &velocityGradAtIP() const
     { return velocityGrad(); }
 
+    /*!
+     * \brief Return the eddy viscosity (if implemented).
+     */
+    const Scalar eddyViscosity() const
+    { return 0; }
+
 //    /*!
 //     * \brief Return the divergence of the normal velocity at the
 //     *        integration point.
diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh
index d53d5c38f3e1947ac4fd65da2e1142ec7d7bfb1f..79929fa8e7156bd687ab1b725b43f5861cfe8aee 100644
--- a/dumux/freeflow/stokes/stokeslocalresidual.hh
+++ b/dumux/freeflow/stokes/stokeslocalresidual.hh
@@ -228,7 +228,7 @@ protected:
             // TODO: dilatation term has to be accounted for in outflow, coupling, neumann
             //            velGradComp[velIdx] += 2./3*fluxVars.velocityDiv;
 
-            velGradComp *= fluxVars.viscosity();
+            velGradComp *= fluxVars.viscosity() + fluxVars.eddyViscosity();
 
             flux[momentumXIdx + velIdx] -=
                 velGradComp*fluxVars.face().normal;
@@ -249,9 +249,8 @@ protected:
         if (calculateNavierStokes)
         {
             for (int dimIndex = 0; dimIndex < dim; ++dimIndex)
-            {
-                flux[momentumXIdx + dimIndex] += up.density() * up.velocity()[dimIndex] * fluxVars.normalVelocity();
-            }
+                flux[momentumXIdx + dimIndex] +=
+                        up.density() * up.velocity()[dimIndex] * fluxVars.normalVelocity();
         }
     }
 
diff --git a/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh b/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh
index f2b9565b9cf8a58871329f124a082f5abbb8a53b..bd1ccef29fcae6fbd1c35789c6780f7f1a944469 100644
--- a/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh
+++ b/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh
@@ -87,7 +87,7 @@ public:
     /*!
      * \brief Return the mass fraction of the transported component at the integration point.
      */
-    Scalar massFraction() const
+    const Scalar massFraction() const
     { return massFraction_; }
 
     /*!
@@ -100,9 +100,15 @@ public:
     /*!
      * \brief Return the molar diffusion coefficient at the integration point.
      */
-    Scalar diffusionCoeff() const
+    const Scalar diffusionCoeff() const
     { return diffusionCoeff_; }
 
+    /*!
+     * \brief Return the eddy diffusivity (if implemented).
+     */
+    const Scalar eddyDiffusivity() const
+    { return 0; }
+
     /*!
      * \brief Return the molar diffusion coefficient at the integration point.
      */
diff --git a/dumux/freeflow/stokes2c/stokes2clocalresidual.hh b/dumux/freeflow/stokes2c/stokes2clocalresidual.hh
index 98e6021c5b28a368427277bad00820a41995c15a..3e934d660c25ff50f925a1efe7267c1670ec8048 100644
--- a/dumux/freeflow/stokes2c/stokes2clocalresidual.hh
+++ b/dumux/freeflow/stokes2c/stokes2clocalresidual.hh
@@ -152,7 +152,7 @@ public:
             flux[transportEqIdx] -=
                 fluxVars.moleFractionGrad()[dimIdx] *
                 fluxVars.face().normal[dimIdx] *
-                fluxVars.diffusionCoeff() *
+                (fluxVars.diffusionCoeff() + fluxVars.eddyDiffusivity()) *
                 fluxVars.molarDensity() *
                 FluidSystem::molarMass(transportCompIdx);
 
diff --git a/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh b/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh
index 66c6a8ce4b09ab39a6938cff7d2c52e5c42cac6a..ba5c448bc2b72165c7a1840aa64fe463305e3ed6 100644
--- a/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh
+++ b/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh
@@ -80,7 +80,7 @@ public:
     /*!
      * \brief Returns the heat conductivity at the integration point.
      */
-    Scalar heatConductivity() const
+    const Scalar heatConductivity() const
     { return heatConductivity_; }
     
     /*!
@@ -96,6 +96,12 @@ public:
     const DimVector &temperatureGrad() const
     { return temperatureGrad_; }
 
+    /*!
+     * \brief Return the eddy conductivity (if implemented).
+     */
+    const Scalar eddyConductivity() const
+    { return 0; }
+
     /*!
      * \brief Returns the temperature gradient at the integration point.
      */
diff --git a/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh b/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh
index e26c0da7cbcbff5e042eabf51abd5589d0c0c6ce..7b21e7a33b7824af7437341e1cc996cd3e34aabb 100644
--- a/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh
+++ b/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh
@@ -137,7 +137,7 @@ public:
             flux[energyEqIdx] -=
                 fluxVars.temperatureGrad()[dimIdx] *
                 fluxVars.face().normal[dimIdx] *
-                fluxVars.heatConductivity();
+                (fluxVars.heatConductivity() + fluxVars.eddyConductivity());
     }
 };