diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
index a080d51cac93a907b8f3c431d90b274d6fb12022..4ce64ab70d272e896d2f1a9a9c6b77aefffefae2 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
@@ -97,11 +97,11 @@ public:
         // calculate advective flux
         auto upwindTermK = [](const auto& volVars)
         {
-            return volVars.turbulentKineticEnergy();
+            return volVars.turbulentKineticEnergy() * volVars.density();
         };
         auto upwindTermEpsilon = [](const auto& volVars)
         {
-            return volVars.dissipationTilde();
+            return volVars.dissipationTilde() * volVars.density();
         };
 
         flux[turbulentKineticEnergyEqIdx]
@@ -116,14 +116,14 @@ public:
         const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
         // effective diffusion coefficients
-        Scalar insideCoeff_k = insideVolVars.kinematicViscosity()
-                               + insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaK();
-        Scalar outsideCoeff_k = outsideVolVars.kinematicViscosity()
-                                + outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaK();
-        Scalar insideCoeff_e = insideVolVars.kinematicViscosity()
-                               + insideVolVars.kinematicEddyViscosity() / insideVolVars.sigmaEpsilon();
-        Scalar outsideCoeff_e = outsideVolVars.kinematicViscosity()
-                                + outsideVolVars.kinematicEddyViscosity() / outsideVolVars.sigmaEpsilon();
+        Scalar insideCoeff_k = insideVolVars.viscosity() + insideVolVars.kinematicEddyViscosity()
+                                                         * insideVolVars.density() / insideVolVars.sigmaK();
+        Scalar outsideCoeff_k = outsideVolVars.viscosity() + outsideVolVars.kinematicEddyViscosity()
+                                                           * outsideVolVars.density() / outsideVolVars.sigmaK();
+        Scalar insideCoeff_e = insideVolVars.viscosity() + insideVolVars.kinematicEddyViscosity()
+                                                         * insideVolVars.density() / insideVolVars.sigmaEpsilon();
+        Scalar outsideCoeff_e = outsideVolVars.viscosity() + outsideVolVars.kinematicEddyViscosity()
+                                                           * outsideVolVars.density() / outsideVolVars.sigmaEpsilon();
 
         // scale by extrusion factor
         insideCoeff_k *= insideVolVars.extrusionFactor();
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh
index b29e5bec5849b7a3ef37a22a9a059eb71a96eadd..1c692aa425236fa03f22e0a32c5dc2b256c2d7ef 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/localresidual.hh
@@ -82,8 +82,8 @@ public:
     {
         CellCenterPrimaryVariables storage = ParentType::computeStorageForCellCenter(problem, scv, volVars);
 
-        storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy();
-        storage[dissipationEqIdx] = volVars.dissipationTilde();
+        storage[turbulentKineticEnergyEqIdx] = volVars.turbulentKineticEnergy() * volVars.density();
+        storage[dissipationEqIdx] = volVars.dissipationTilde() * volVars.density();
 
         return storage;
     }
@@ -101,22 +101,22 @@ public:
         const auto& volVars = elemVolVars[scv];
 
         // production
-        source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.kinematicEddyViscosity()
+        source[turbulentKineticEnergyEqIdx] += 2.0 * volVars.dynamicEddyViscosity()
                                                * volVars.stressTensorScalarProduct();
         source[dissipationEqIdx] += volVars.cOneEpsilon() * volVars.fOne()
                                     * volVars.dissipationTilde() / volVars.turbulentKineticEnergy()
-                                    * 2.0 * volVars.kinematicEddyViscosity()
+                                    * 2.0 * volVars.dynamicEddyViscosity()
                                     * volVars.stressTensorScalarProduct();
 
         // destruction
-        source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde();
-        source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo()
+        source[turbulentKineticEnergyEqIdx] -= volVars.dissipationTilde() * volVars.density();
+        source[dissipationEqIdx] -= volVars.cTwoEpsilon() * volVars.fTwo() * volVars.density()
                                     * volVars.dissipationTilde() * volVars.dissipationTilde()
                                     / volVars.turbulentKineticEnergy();
 
         // dampening functions
-        source[turbulentKineticEnergyEqIdx] -= volVars.dValue();
-        source[dissipationEqIdx] += volVars.eValue();
+        source[turbulentKineticEnergyEqIdx] -= volVars.dValue() * volVars.density();
+        source[dissipationEqIdx] += volVars.eValue() * volVars.density();
 
         return source;
     }