diff --git a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index 18f8d49995241dda50bf28a0bde304a2335b2e9f..f1dc4e9203767f045707e99ffb66103d7e5e47e4 100644
--- a/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -168,10 +168,10 @@ LaunderSharmaKE::LaunderSharmaKE
             "mut",
             runTime_.timeName(),
             mesh_,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        rho_*Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_)
+        autoCreateLowReMut("mut", mesh_)
     ),
 
     alphat_
@@ -187,6 +187,9 @@ LaunderSharmaKE::LaunderSharmaKE
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    mut_ = rho_*Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
+
     alphat_ = mut_/Prt_;
     alphat_.correctBoundaryConditions();
 
diff --git a/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C b/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
index 842d9125ea27d823cfcfa644dd0a4f0c0f775d22..d2e291d915de4c2c31a6b3e4ac84797827dacb13 100644
--- a/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
+++ b/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
@@ -29,6 +29,7 @@ License
 #include "calculatedFvPatchField.H"
 #include "alphatWallFunctionFvPatchScalarField.H"
 #include "mutkWallFunctionFvPatchScalarField.H"
+#include "mutLowReWallFunctionFvPatchScalarField.H"
 #include "epsilonWallFunctionFvPatchScalarField.H"
 #include "kqRWallFunctionFvPatchField.H"
 #include "omegaWallFunctionFvPatchScalarField.H"
@@ -182,6 +183,76 @@ tmp<volScalarField> autoCreateMut
 }
 
 
+tmp<volScalarField> autoCreateLowReMut
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    IOobject mutHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (mutHeader.headerOk())
+    {
+        return tmp<volScalarField>(new volScalarField(mutHeader, mesh));
+    }
+    else
+    {
+        Info<< "--> Creating " << fieldName
+            << " to employ run-time selectable wall functions" << endl;
+
+        const fvBoundaryMesh& bm = mesh.boundary();
+
+        wordList mutBoundaryTypes(bm.size());
+
+        forAll(bm, patchI)
+        {
+            if (isA<wallFvPatch>(bm[patchI]))
+            {
+                mutBoundaryTypes[patchI] =
+                    RASModels::mutLowReWallFunctionFvPatchScalarField::typeName;
+            }
+            else
+            {
+                mutBoundaryTypes[patchI] =
+                    calculatedFvPatchField<scalar>::typeName;
+            }
+        }
+
+        tmp<volScalarField> mut
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh,
+                dimensionedScalar("zero", dimDensity*dimArea/dimTime, 0.0),
+                mutBoundaryTypes
+            )
+        );
+
+        Info<< "    Writing new " << fieldName << endl;
+        mut().write();
+
+        return mut;
+    }
+}
+
+
 tmp<volScalarField> autoCreateEpsilon
 (
     const word& fieldName,
diff --git a/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H b/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H
index d36be7583cc2f069f00a7b95aaccc14b48a20a4e..31e38aa4152767481aad7a53351e479c7080c51b 100644
--- a/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H
+++ b/src/turbulenceModels/compressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H
@@ -53,6 +53,13 @@ namespace compressible
         const fvMesh& mesh
     );
 
+    //- mut for Low-Reynolds number models
+    tmp<volScalarField> autoCreateLowReMut
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
     //- alphat
     tmp<volScalarField> autoCreateAlphat
     (
diff --git a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
index 321a7b5ff91df01571d207d804f8f2ed165bfde7..1edcb224b556829a6a0aa2a1a8ce383ed251b662 100644
--- a/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
@@ -127,8 +127,22 @@ LamBremhorstKE::LamBremhorstKE
        *(scalar(1) + 20.5/(Rt_ + SMALL))
     ),
 
-    nut_(Cmu_*fMu_*sqr(k_)/(epsilon_ + epsilonSmall_))
+    nut_
+    (
+        IOobject
+        (
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateLowReNut("nut", mesh_)
+    )
 {
+    nut_ = Cmu_*fMu_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    nut_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
diff --git a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
index ad2ebe0821eb994026b69df4fded8daf31fa70d5..3bbdcb6a4734a47eee3a0543d5ed6d741d653e12 100644
--- a/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
+++ b/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
@@ -133,8 +133,22 @@ LaunderSharmaKE::LaunderSharmaKE
         mesh_
     ),
 
-    nut_(Cmu_*fMu()*sqr(k_)/(epsilonTilda_ + epsilonSmall_))
+    nut_
+    (
+        IOobject
+        (
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateLowReNut("nut", mesh_)
+    )
 {
+    nut_ = Cmu_*fMu()*sqr(k_)/(epsilonTilda_ + epsilonSmall_);
+    nut_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
index 2fdb30e8cf9fbfe68c7e602ee8ad566988b53234..2146553f0f6b758230156e68243fdb28af9b598c 100644
--- a/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C
@@ -28,6 +28,8 @@ License
 #include "wallFvPatch.H"
 #include "addToRunTimeSelectionTable.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -233,20 +235,16 @@ LienCubicKELowRe::LienCubicKELowRe
 
     nut_
     (
-        Cmu_
-       *(
-            scalar(1) - exp(-Am_*yStar_))
-           /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL
-        )
-       *sqr(k_)/(epsilon_ + epsilonSmall_)
-        // cubic term C5, implicit part
-      + max
+        IOobject
         (
-            C5viscosity_,
-            dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
-        )
+            "nut",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateLowReNut("nut", mesh_)
     ),
-    // turbulent viscosity, with implicit part of C5
 
     nonlinearStress_
     (
@@ -282,6 +280,21 @@ LienCubicKELowRe::LienCubicKELowRe
         )
     )
 {
+    nut_ = Cmu_
+       *(
+            scalar(1) - exp(-Am_*yStar_))
+           /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL
+        )
+       *sqr(k_)/(epsilon_ + epsilonSmall_)
+        // cubic term C5, implicit part
+      + max
+        (
+            C5viscosity_,
+            dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
+        );
+
+    nut_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
index b7a82d6d1bdd788b3ce5bb598e677eae27560f10..6e47e1d926f77bd2e5a794decc8a6c4c16ddfd65 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@@ -28,6 +28,8 @@ License
 #include "wallFvPatch.H"
 #include "addToRunTimeSelectionTable.H"
 
+#include "backwardsCompatibilityWallFunctions.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -167,11 +169,22 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
 
     nut_
     (
-        Cmu_*(scalar(1) - exp(-Am_*yStar_))
-       /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL)*sqr(k_)
-       /(epsilon_ + epsilonSmall_)
+        IOobject
+        (
+            "epsilon",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateLowReNut("nut", mesh_)
     )
 {
+    nut_ = Cmu_*(scalar(1) - exp(-Am_*yStar_))
+       /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL)*sqr(k_)
+       /(epsilon_ + epsilonSmall_);
+    nut_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
diff --git a/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C b/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
index b2bcedff458b7d5311dc6966ba923d18d475cc3a..3de77712fdcfa64f1aa6e8a90b37167e2e902583 100644
--- a/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
+++ b/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
@@ -28,6 +28,7 @@ License
 
 #include "calculatedFvPatchField.H"
 #include "nutkWallFunctionFvPatchScalarField.H"
+#include "nutLowReWallFunctionFvPatchScalarField.H"
 #include "epsilonWallFunctionFvPatchScalarField.H"
 #include "kqRWallFunctionFvPatchField.H"
 #include "omegaWallFunctionFvPatchScalarField.H"
@@ -111,6 +112,76 @@ tmp<volScalarField> autoCreateNut
 }
 
 
+tmp<volScalarField> autoCreateLowReNut
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    IOobject nutHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (nutHeader.headerOk())
+    {
+        return tmp<volScalarField>(new volScalarField(nutHeader, mesh));
+    }
+    else
+    {
+        Info<< "--> Creating " << fieldName
+            << " to employ run-time selectable wall functions" << endl;
+
+        const fvBoundaryMesh& bm = mesh.boundary();
+
+        wordList nutBoundaryTypes(bm.size());
+
+        forAll(bm, patchI)
+        {
+            if (isA<wallFvPatch>(bm[patchI]))
+            {
+                nutBoundaryTypes[patchI] =
+                    RASModels::nutLowReWallFunctionFvPatchScalarField::typeName;
+            }
+            else
+            {
+                nutBoundaryTypes[patchI] =
+                    calculatedFvPatchField<scalar>::typeName;
+            }
+        }
+
+        tmp<volScalarField> nut
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh,
+                dimensionedScalar("zero", dimArea/dimTime, 0.0),
+                nutBoundaryTypes
+            )
+        );
+
+        Info<< "    Writing new " << fieldName << endl;
+        nut().write();
+
+        return nut;
+    }
+}
+
+
 tmp<volScalarField> autoCreateEpsilon
 (
     const word& fieldName,
diff --git a/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H b/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H
index 615834f9e804b54b291db21de1dae1269e0d7ada..e3d2c738851b3cb153c3b5e55382815134479955 100644
--- a/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H
+++ b/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.H
@@ -53,6 +53,13 @@ namespace incompressible
         const fvMesh& mesh
     );
 
+    //- nut for Low-Reynolds number models
+    tmp<volScalarField> autoCreateLowReNut
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
     //- epsilon
     tmp<volScalarField> autoCreateEpsilon
     (