diff --git a/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
index 85ff9c44..57a6bca6 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
@@ -93,14 +93,14 @@ void epsilonLowReWallFunctionFvPatchScalarField::calculate
 
         if (yPlus > yPlusLam_)
         {
-            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+            epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
         }
         else
         {
-            epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
+            epsilon[cellI] += w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
         }
 
-        G[cellI] =
+        G[cellI] +=
             w
            *(nutw[faceI] + nuw[faceI])
            *magGradUw[faceI]
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
index 347993c1..1cab0ae6 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
@@ -98,15 +98,15 @@ void epsilonLowReWallFunctionFvPatchScalarField::calculate
 
         if (yPlus > yPlusLam_)
         {
-            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+            epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
         }
         else
         {
-            epsilon[cellI] =
+            epsilon[cellI] +=
                 w*2.0*k[cellI]*muw[faceI]/rhow[faceI]/sqr(y[faceI]);
         }
 
-        G[cellI] =
+        G[cellI] +=
             w
            *(mutw[faceI] + muw[faceI])
            *magGradUw[faceI]
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
index d9898bed..93946739 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
@@ -96,14 +96,14 @@ void epsilonLowReWallFunctionFvPatchScalarField::calculate
 
         if (yPlus > yPlusLam_)
         {
-            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
+            epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
         }
         else
         {
-            epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
+            epsilon[cellI] += w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
         }
 
-        G[cellI] =
+        G[cellI] +=
             w
            *(nutw[faceI] + nuw[faceI])
            *magGradUw[faceI]