diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
index 9a721c9f1ffe02747df166bea6fb09b31074cbd8..71a3cb01c73b95ac549323114283cbe3f3ffff53 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -238,21 +238,24 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
 
         const scalar w = cornerWeights[facei];
 
-        if (yPlus > yPlusLam_)
+        // Default high-Re form
+        scalar epsilonc = w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
+        scalar Gc =
+            w
+           *(nutw[facei] + nuw[facei])
+           *magGradUw[facei]
+           *Cmu25*sqrt(k[celli])
+           /(kappa_*y[facei]);
+
+        if (lowReCorrection_ && yPlus < yPlusLam_)
         {
-            epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
-
-            G0[celli] +=
-                w
-               *(nutw[facei] + nuw[facei])
-               *magGradUw[facei]
-               *Cmu25*sqrt(k[celli])
-               /(kappa_*y[facei]);
-        }
-        else
-        {
-            epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
+            epsilonc = w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
+            Gc = 0;
         }
+
+        epsilon0[celli] += epsilonc;
+
+        G0[celli] += Gc;
     }
 }
 
@@ -273,6 +276,7 @@ epsilonWallFunctionFvPatchScalarField
     yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
     G_(),
     epsilon_(),
+    lowReCorrection_(false),
     initialised_(false),
     master_(-1),
     cornerWeights_()
@@ -297,6 +301,7 @@ epsilonWallFunctionFvPatchScalarField
     yPlusLam_(ptf.yPlusLam_),
     G_(),
     epsilon_(),
+    lowReCorrection_(ptf.lowReCorrection_),
     initialised_(false),
     master_(-1),
     cornerWeights_()
@@ -320,6 +325,7 @@ epsilonWallFunctionFvPatchScalarField
     yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
     G_(),
     epsilon_(),
+    lowReCorrection_(dict.lookupOrDefault("lowReCorrection", false)),
     initialised_(false),
     master_(-1),
     cornerWeights_()
@@ -344,6 +350,7 @@ epsilonWallFunctionFvPatchScalarField
     yPlusLam_(ewfpsf.yPlusLam_),
     G_(),
     epsilon_(),
+    lowReCorrection_(ewfpsf.lowReCorrection_),
     initialised_(false),
     master_(-1),
     cornerWeights_()
@@ -366,6 +373,7 @@ epsilonWallFunctionFvPatchScalarField
     yPlusLam_(ewfpsf.yPlusLam_),
     G_(),
     epsilon_(),
+    lowReCorrection_(ewfpsf.lowReCorrection_),
     initialised_(false),
     master_(-1),
     cornerWeights_()
@@ -492,11 +500,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::updateWeightedCoeffs
 
     typedef DimensionedField<scalar, volMesh> FieldType;
 
-    FieldType& G =
-        const_cast<FieldType&>
-        (
-            db().lookupObject<FieldType>(turbModel.GName())
-        );
+    FieldType& G = db().lookupObjectRef<FieldType>(turbModel.GName());
 
     FieldType& epsilon = const_cast<FieldType&>(internalField());
 
@@ -560,7 +564,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::manipulateMatrix
 
     forAll(weights, facei)
     {
-        // Anly set the values if the weights are > tolerance
+        // Only set the values if the weights are > tolerance
         if (weights[facei] > tolerance_)
         {
             nConstrainedCells++;
@@ -594,6 +598,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     writeLocalEntries(os);
     fixedValueFvPatchField<scalar>::write(os);
+    os.writeEntry("lowReCorrection", lowReCorrection_);
 }
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
index f75aa95b2379770e01bc6f010c8f34d464c22b98..b52bd48aa120babb21a50733be7bedc66e6ff079 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
+++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H
@@ -42,15 +42,19 @@ Description
         G       | turblence generation field
     \endvartable
 
-    The model switches between laminar and turbulent functions based on the
-    laminar-to-turbulent y+ value derived from kappa and E.
+    The low-Re correction is activated by setting the entry
+    \c lowReCorrection to 'on'; in this mode the model switches between
+    laminar and turbulent functions based on the laminar-to-turbulent y+ value
+    derived from kappa and E.  When the \c lowReCorrection is inactive, the
+    wall function operates in high-Re mode.
 
 Usage
     \table
-        Property     | Description             | Required    | Default value
-        Cmu          | model coefficient       | no          | 0.09
-        kappa        | Von Karman constant     | no          | 0.41
-        E            | model coefficient       | no          | 9.8
+        Property        | Description              | Required  | Default value
+        Cmu             | model coefficient        | no        | 0.09
+        kappa           | Von Karman constant      | no        | 0.41
+        E               | model coefficient        | no        | 9.8
+        lowReCorrection | Low-Re correction active | no        | off
     \endtable
 
     Example of the boundary condition specification:
@@ -115,6 +119,9 @@ protected:
         //- Local copy of turbulence epsilon field
         scalarField epsilon_;
 
+        //- Apply low-Re correction term; default = no
+        bool lowReCorrection_;
+
         //- Initialised flag
         bool initialised_;