diff --git a/src/finiteVolume/cfdTools/general/bound/bound.C b/src/finiteVolume/cfdTools/general/bound/bound.C
index 7054434c0f..93dac0faf9 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.C
+++ b/src/finiteVolume/cfdTools/general/bound/bound.C
@@ -30,7 +30,8 @@ License
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-void Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
+Foam::volScalarField&
+Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
 {
     const scalar minVsf = min(vsf).value();
 
@@ -55,6 +56,8 @@ void Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
 
         vsf.boundaryField() = max(vsf.boundaryField(), lowerBound.value());
     }
+
+    return vsf;
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/bound/bound.H b/src/finiteVolume/cfdTools/general/bound/bound.H
index f301c5a217..1f6418ce5c 100644
--- a/src/finiteVolume/cfdTools/general/bound/bound.H
+++ b/src/finiteVolume/cfdTools/general/bound/bound.H
@@ -50,8 +50,9 @@ namespace Foam
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
 //- Bound the given scalar field if it has gone unbounded.
+//  Return the bounded field.
 //  Used extensively in RAS and LES turbulence models.
-void bound(volScalarField&, const dimensionedScalar& lowerBound);
+volScalarField& bound(volScalarField&, const dimensionedScalar& lowerBound);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
index 5f036cff10..5212e7bbad 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
@@ -119,7 +119,7 @@ LamBremhorstKE::LamBremhorstKE
 
     y_(mesh_),
 
-    Rt_(sqr(k_)/(nu()*epsilon_)),
+    Rt_(sqr(k_)/(nu()*bound(epsilon_, epsilonMin_))),
 
     fMu_
     (
@@ -141,7 +141,7 @@ LamBremhorstKE::LamBremhorstKE
     )
 {
     bound(k_, kMin_);
-    bound(epsilon_, epsilonMin_);
+    // already bounded: bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*fMu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
index d91f974fbe..c5c783f125 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
@@ -173,12 +173,12 @@ LienCubicKE::LienCubicKE
     gradU_(fvc::grad(U)),
     eta_
     (
-        k_/(epsilon_ + epsilonMin_)
+        k_/bound(epsilon_, epsilonMin_)
        *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
     ),
     ksi_
     (
-        k_/(epsilon_ + epsilonMin_)
+        k_/epsilon_
        *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
     ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
@@ -235,7 +235,7 @@ LienCubicKE::LienCubicKE
     )
 {
     bound(k_, kMin_);
-    bound(epsilon_, epsilonMin_);
+    // already bounded: bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;
     nut_.correctBoundaryConditions();
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
index 046f5a6fa1..afd85ff475 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
@@ -222,12 +222,12 @@ LienCubicKELowRe::LienCubicKELowRe
     gradU_(fvc::grad(U)),
     eta_
     (
-        k_/(epsilon_ + epsilonMin_)
+        k_/bound(epsilon_, epsilonMin_)
        *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
     ),
     ksi_
     (
-        k_/(epsilon_ + epsilonMin_)
+        k_/epsilon_
        *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
     ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
@@ -235,7 +235,7 @@ LienCubicKELowRe::LienCubicKELowRe
 
     C5viscosity_
     (
-        -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_ + epsilonMin_)
+        -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()))
     ),
 
@@ -260,7 +260,7 @@ LienCubicKELowRe::LienCubicKELowRe
         symm
         (
             // quadratic terms
-            pow3(k_)/sqr(epsilon_ + epsilonMin_)
+            pow3(k_)/sqr(epsilon_)
            *(
                 Ctau1_/fEta_
                *(
@@ -271,7 +271,7 @@ LienCubicKELowRe::LienCubicKELowRe
               + Ctau3_/fEta_*(gradU_.T() & gradU_)
             )
             // cubic term C4
-          - 20.0*pow4(k_)/pow3(epsilon_ + epsilonMin_)
+          - 20.0*pow4(k_)/pow3(epsilon_)
            *pow3(Cmu_)
            *(
                 ((gradU_ & gradU_) & gradU_.T())
@@ -289,7 +289,7 @@ LienCubicKELowRe::LienCubicKELowRe
     )
 {
     bound(k_, kMin_);
-    bound(epsilon_, epsilonMin_);
+    // already bounded: bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_
       * (scalar(1) - exp(-Am_*yStar_))
diff --git a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
index c885d87f16..fdcae3d187 100644
--- a/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
+++ b/src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
@@ -192,25 +192,25 @@ NonlinearKEShih::NonlinearKEShih
     gradU_(fvc::grad(U)),
     eta_
     (
-        k_/(epsilon_ + epsilonMin_)
+        k_/bound(epsilon_, epsilonMin_)
        *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
     ),
     ksi_
     (
-        k_/(epsilon_+ epsilonMin_)
+        k_/epsilon_
        *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
     ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
     fEta_(A2_ + pow(eta_, 3.0)),
 
-    nut_("nut", Cmu_*sqr(k_)/(epsilon_ + epsilonMin_)),
+    nut_("nut", Cmu_*sqr(k_)/epsilon_),
 
     nonlinearStress_
     (
         "nonlinearStress",
         symm
         (
-            pow3(k_)/sqr(epsilon_ + epsilonMin_)
+            pow3(k_)/sqr(epsilon_)
            *(
                 Ctau1_/fEta_
                *(
@@ -224,7 +224,7 @@ NonlinearKEShih::NonlinearKEShih
     )
 {
     bound(k_, kMin_);
-    bound(epsilon_, epsilonMin_);
+    // already bounded: bound(epsilon_, epsilonMin_);
 
     #include "wallNonlinearViscosityI.H"
 
diff --git a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
index 422f8b1463..41b1114429 100644
--- a/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
+++ b/src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
@@ -179,7 +179,7 @@ qZeta::qZeta
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        epsilon_/(2.0*q_),
+        epsilon_/(2.0*bound(q_, qMin_)),
         epsilon_.boundaryField().types()
     ),
 
@@ -198,7 +198,7 @@ qZeta::qZeta
 {
     bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
-    bound(q_, qMin_);
+    // already bounded: bound(q_, qMin_);
     bound(zeta_, zetaMin_);
 
     nut_ = Cmu_*fMu()*sqr(k_)/epsilon_;