diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
index afe97f7f7f38c406a67fd0d66c9c42e38735e783..4a2a35635bbc3ffe23989b475322249d53617428 100644
--- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
+++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
@@ -130,7 +130,7 @@ void PDRkEpsilon::correct()
     tgradU.clear();
 
     // Update espsilon and G at the wall
-    epsilon_.boundaryField().updateCoeffs();
+    epsilon_.boundaryFieldRef().updateCoeffs();
 
     // Add the blockage generation term so that it is included consistently
     // in both the k and epsilon equations
@@ -163,7 +163,7 @@ void PDRkEpsilon::correct()
 
     epsEqn.ref().relax();
 
-    epsEqn.ref().boundaryManipulate(epsilon_.boundaryField());
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
 
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);
diff --git a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
index 0208818ee0034cb202e0a4e1a5d213c8c2a17972..6ca24555b4139228603ff94d73aef5f07de68ec6 100644
--- a/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
+++ b/applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
@@ -119,9 +119,11 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
         }
     }
 
+    volScalarField::GeometricBoundaryField& xieqBf = xieq.boundaryFieldRef();
+
     forAll(xieq.boundaryField(), patchi)
     {
-        scalarField& xieqp = xieq.boundaryField()[patchi];
+        scalarField& xieqp = xieqBf[patchi];
         const scalarField& Kp = K.boundaryField()[patchi];
         const scalarField& Map = Ma.boundaryField()[patchi];
         const scalarField& upBySup = upBySu.boundaryField()[patchi];
diff --git a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
index fd07b0b3cbfcada4c87eda806251c708cdaf73b8..17f92a85f214d32cbfd6bf9d8ce561f4bd3cfcea 100644
--- a/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
+++ b/applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
@@ -266,9 +266,11 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi);
     }
 
-    forAll(Su0.boundaryField(), patchi)
+    volScalarField::GeometricBoundaryField& Su0Bf = Su0.boundaryFieldRef();
+
+    forAll(Su0Bf, patchi)
     {
-        scalarField& Su0p = Su0.boundaryField()[patchi];
+        scalarField& Su0p = Su0Bf[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
 
@@ -313,9 +315,11 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi[celli]);
     }
 
-    forAll(Su0.boundaryField(), patchi)
+    volScalarField::GeometricBoundaryField& Su0Bf = Su0.boundaryFieldRef();
+
+    forAll(Su0Bf, patchi)
     {
-        scalarField& Su0p = Su0.boundaryField()[patchi];
+        scalarField& Su0p = Su0Bf[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
@@ -365,9 +369,11 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
         ma[celli] = Ma(phi[celli]);
     }
 
-    forAll(ma.boundaryField(), patchi)
+    volScalarField::GeometricBoundaryField& maBf = ma.boundaryFieldRef();
+
+    forAll(maBf, patchi)
     {
-        scalarField& map = ma.boundaryField()[patchi];
+        scalarField& map = maBf[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
 
         forAll(map, facei)
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
index 3b43b5dcec39083ca06b08ecc07871b5c418924a..fb0fbc7079565c6a47d65d3c76117c6bc6b91a8f 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
@@ -189,7 +189,7 @@ int main(int argc, char *argv[])
             rhoU.dimensionedInternalField()
            /rho.dimensionedInternalField();
         U.correctBoundaryConditions();
-        rhoU.boundaryField() == rho.boundaryField()*U.boundaryField();
+        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
 
         if (!inviscid)
         {
@@ -223,7 +223,7 @@ int main(int argc, char *argv[])
         e = rhoE/rho - 0.5*magSqr(U);
         e.correctBoundaryConditions();
         thermo.correct();
-        rhoE.boundaryField() ==
+        rhoE.boundaryFieldRef() ==
             rho.boundaryField()*
             (
                 e.boundaryField() + 0.5*magSqr(U.boundaryField())
@@ -244,7 +244,7 @@ int main(int argc, char *argv[])
             rho.dimensionedInternalField()
            /psi.dimensionedInternalField();
         p.correctBoundaryConditions();
-        rho.boundaryField() == psi.boundaryField()*p.boundaryField();
+        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
 
         turbulence->correct();
 
diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
index 3754632d3d4be9724dbc7544c96497c2acf68889..083a1b1b5c1ebb2ed0185e436ea1559a6a295838 100644
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
@@ -182,7 +182,7 @@ int main(int argc, char *argv[])
             rhoU.dimensionedInternalField()
            /rho.dimensionedInternalField();
         U.correctBoundaryConditions();
-        rhoU.boundaryField() == rho.boundaryField()*U.boundaryField();
+        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();
 
         if (!inviscid)
         {
@@ -216,7 +216,7 @@ int main(int argc, char *argv[])
         e = rhoE/rho - 0.5*magSqr(U);
         e.correctBoundaryConditions();
         thermo.correct();
-        rhoE.boundaryField() ==
+        rhoE.boundaryFieldRef() ==
             rho.boundaryField()*
             (
                 e.boundaryField() + 0.5*magSqr(U.boundaryField())
@@ -237,7 +237,7 @@ int main(int argc, char *argv[])
             rho.dimensionedInternalField()
            /psi.dimensionedInternalField();
         p.correctBoundaryConditions();
-        rho.boundaryField() == psi.boundaryField()*p.boundaryField();
+        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();
 
         turbulence->correct();
 
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
index 14bbbe8708650e3aaab85a355e759832387ce386..db4b5972e98562193e493dbea7c28fff5b368851 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C
@@ -923,7 +923,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
 
-    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryField());
+    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index 9d6dd2dc18b68f8e2adbd1b422ad2a54c997d915..f75f4d0f51e6d8f1aafaf5e58e351fbec488c860 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -54,11 +54,14 @@
         phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
     }
 
+    surfaceScalarField::GeometricBoundaryField& phicBf =
+        phic.boundaryFieldRef();
+
     // Do not compress interface at non-coupled boundary faces
     // (inlets, outlets etc.)
     forAll(phic.boundaryField(), patchi)
     {
-        fvsPatchScalarField& phicp = phic.boundaryField()[patchi];
+        fvsPatchScalarField& phicp = phicBf[patchi];
 
         if (!phicp.coupled())
         {
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
index 8343ceb79038e8d3342b3746f97e9546e366ebcc..d5ee2f25f872d06aca73833b6ba2dd7d43bda215 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -134,7 +134,7 @@ void Foam::threePhaseInterfaceProperties::calculateK()
     // Face unit interface normal
     surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
 
-    correctContactAngle(nHatfv.boundaryField());
+    correctContactAngle(nHatfv.boundaryFieldRef());
 
     // Face unit interface normal flux
     nHatf_ = nHatfv & Sf;
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 3f514904eca92299229768cf0e69cba88a90fa1a..62e2a97a3016077f74a430e335ad119009874976 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -124,11 +124,13 @@ void Foam::multiphaseSystem::solveAlphas()
             );
         }
 
+        surfaceScalarField::GeometricBoundaryField& alphaPhiCorrBf =
+            alphaPhiCorr.boundaryFieldRef();
+
         // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhiCorr.boundaryField(), patchi)
+        forAll(alphaPhiCorrBf, patchi)
         {
-            fvsPatchScalarField& alphaPhiCorrp =
-                alphaPhiCorr.boundaryField()[patchi];
+            fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];
 
             if (!alphaPhiCorrp.coupled())
             {
@@ -372,7 +374,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
 
-    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryField());
+    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
@@ -666,6 +668,9 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
         }
     }
 
+    volVectorField::GeometricBoundaryField& SvmBf =
+        tSvm.ref().boundaryFieldRef();
+
     // Remove virtual mass at fixed-flux boundaries
     forAll(phase.phi().boundaryField(), patchi)
     {
@@ -677,7 +682,7 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
             )
         )
         {
-            tSvm.ref().boundaryField()[patchi] = Zero;
+            SvmBf[patchi] = Zero;
         }
     }
 
@@ -713,6 +718,8 @@ Foam::multiphaseSystem::dragCoeffs() const
                 )
             ).ptr();
 
+        volScalarField::GeometricBoundaryField& Kbf = Kptr->boundaryFieldRef();
+
         // Remove drag at fixed-flux boundaries
         forAll(dm.phase1().phi().boundaryField(), patchi)
         {
@@ -724,7 +731,7 @@ Foam::multiphaseSystem::dragCoeffs() const
                 )
             )
             {
-                Kptr->boundaryField()[patchi] = 0.0;
+                Kbf[patchi] = 0.0;
             }
         }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index bbcc36f4e772114c866bd51cd5210956e51bf5b7..f65fdf8de168eb17070f3487cbb2d8ab2349a8fe 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -203,7 +203,7 @@
 
         setSnGrad<fixedFluxPressureFvPatchScalarField>
         (
-            p_rgh.boundaryField(),
+            p_rgh.boundaryFieldRef(),
             (
                 phiHbyA.boundaryField() - MRF.relative(phib)
             )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index ea08ba29865801976b32ffb25801466fad2b786c..f1fed0103d1bbf3a40dbdf5ea16012ee53e642d8 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -523,7 +523,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
 
-    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryField());
+    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C
index a2fbcd0bbcc61d9068ccfbd42e33fbd5c59e82af..886da7c20094725fc741d1157bebce9f2630e68d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/saturationModels/polynomial/polynomial.C
@@ -118,9 +118,11 @@ Foam::saturationModels::polynomial::Tsat
         Tsat[celli] = C_.value(p[celli]);
     }
 
+    volScalarField::GeometricBoundaryField& TsatBf = Tsat.boundaryFieldRef();
+
     forAll(Tsat.boundaryField(), patchi)
     {
-        scalarField& Tsatp = Tsat.boundaryField()[patchi];
+        scalarField& Tsatp = TsatBf[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
 
         forAll(Tsatp, facei)
diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
index a2a95cce05733d902c263dac63511a67309d99d5..784f323e0ca9cf3cd51f29c7d95055fc5c620e45 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C
@@ -50,11 +50,13 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModel::zeroGradWalls
     volVectorField& Fi = tFi.ref();
     const fvPatchList& patches =  Fi.mesh().boundary();
 
+    volVectorField::GeometricBoundaryField& FiBf = Fi.boundaryFieldRef();
+
     forAll(patches, patchi)
     {
         if (isA<wallFvPatch>(patches[patchi]))
         {
-            fvPatchVectorField& Fiw = Fi.boundaryField()[patchi];
+            fvPatchVectorField& Fiw = FiBf[patchi];
             Fiw = Fiw.patchInternalField();
         }
     }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
index f216859b08e0938ebf566e1bae06d68d94a27dd9..a19b995bc598920b08dd044d3badffa911ace8b8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -36,6 +36,9 @@ void Foam::BlendedInterfacialModel<ModelType>::correctFixedFluxBCs
     GeometricField& field
 ) const
 {
+    typename GeometricField::GeometricBoundaryField& fieldBf =
+        field.boundaryFieldRef();
+
     forAll(phase1_.phi()().boundaryField(), patchi)
     {
         if
@@ -46,7 +49,7 @@ void Foam::BlendedInterfacialModel<ModelType>::correctFixedFluxBCs
             )
         )
         {
-            field.boundaryField()[patchi] = Zero;
+            fieldBf[patchi] = Zero;
         }
     }
 }
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index ff1e0b0998e785a62760e51db2c1fe92d0d56fde..8ba5a4a502486e5b3dfb011837b6035b144d304f 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -353,9 +353,9 @@ void Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctThermo()
 
         // Limit the H[12] boundary field to avoid /0
         const scalar HLimit = 1e-4;
-        H1.boundaryField() =
+        H1.boundaryFieldRef() =
             max(H1.boundaryField(), phase1.boundaryField()*HLimit);
-        H2.boundaryField() =
+        H2.boundaryFieldRef() =
             max(H2.boundaryField(), phase2.boundaryField()*HLimit);
 
         volScalarField mDotL
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index b81e174db878bfe6fbd061764179cc07e0c7ed19..6504f943a0e482f41ba84a4487ddba960eb261c5 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -134,11 +134,13 @@ void Foam::multiphaseSystem::solveAlphas()
             );
         }
 
+        surfaceScalarField::GeometricBoundaryField& alphaPhiCorrBf =
+            alphaPhiCorr.boundaryFieldRef();
+
         // Ensure that the flux at inflow BCs is preserved
         forAll(alphaPhiCorr.boundaryField(), patchi)
         {
-            fvsPatchScalarField& alphaPhiCorrp =
-                alphaPhiCorr.boundaryField()[patchi];
+            fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];
 
             if (!alphaPhiCorrp.coupled())
             {
@@ -477,7 +479,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
 {
     tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
 
-    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryField());
+    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
 
     // Simple expression for curvature
     return -fvc::div(tnHatfv & mesh_.Sf());
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 05a049b4c6e2aa4b4e3993ce812c108f2aad68f0..e078414a20b52ef57e912e88e48f32c67ce03ecd 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -182,6 +182,9 @@ while (pimple.correct())
         );
         surfaceScalarField phiCorrCoeff(pos(alphafBar - 0.99));
 
+        surfaceScalarField::GeometricBoundaryField& phiCorrCoeffBf =
+            phiCorrCoeff.boundaryFieldRef();
+
         forAll(mesh.boundary(), patchi)
         {
             // Set ddtPhiCorr to 0 on non-coupled boundaries
@@ -191,7 +194,7 @@ while (pimple.correct())
              || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
             )
             {
-                phiCorrCoeff.boundaryField()[patchi] = 0;
+                phiCorrCoeffBf[patchi] = 0;
             }
         }
 
@@ -281,7 +284,7 @@ while (pimple.correct())
 
         setSnGrad<fixedFluxPressureFvPatchScalarField>
         (
-            p_rgh.boundaryField(),
+            p_rgh.boundaryFieldRef(),
             (
                 phiHbyA.boundaryField() - phib
             )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
index ba312ab248e5c8401ba3311a10caa303a3e74758..008134c21a7cf1f225d3d79a83c3d26c6902a4a9 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H
@@ -149,17 +149,25 @@ while (pimple.correct())
     surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99));
     surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar));
 
-    forAll(mesh.boundary(), patchi)
     {
-        // Set ddtPhiCorr to 0 on non-coupled boundaries
-        if
-        (
-            !mesh.boundary()[patchi].coupled()
-         || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
-        )
+        surfaceScalarField::GeometricBoundaryField& phiCorrCoeff1Bf =
+            phiCorrCoeff1.boundaryFieldRef();
+
+        surfaceScalarField::GeometricBoundaryField& phiCorrCoeff2Bf =
+            phiCorrCoeff2.boundaryFieldRef();
+
+        forAll(mesh.boundary(), patchi)
         {
-            phiCorrCoeff1.boundaryField()[patchi] = 0;
-            phiCorrCoeff2.boundaryField()[patchi] = 0;
+            // Set ddtPhiCorr to 0 on non-coupled boundaries
+            if
+            (
+               !mesh.boundary()[patchi].coupled()
+             || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
+            )
+            {
+                phiCorrCoeff1Bf[patchi] = 0;
+                phiCorrCoeff2Bf[patchi] = 0;
+            }
         }
     }
 
@@ -215,7 +223,7 @@ while (pimple.correct())
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
     (
-        p_rgh.boundaryField(),
+        p_rgh.boundaryFieldRef(),
         (
             phiHbyA.boundaryField()
           - (
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
index 444f2de81f81692ceda80ebb1a0b30f8244c0a17..01c120d481f513aed2fb58ac82cbd7179339ef84 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H
@@ -201,7 +201,7 @@ while (pimple.correct())
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
     (
-        p_rgh.boundaryField(),
+        p_rgh.boundaryFieldRef(),
         (
             phiHbyA.boundaryField()
           - (
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
index 91b49b30e0ea6cf577e50ddb0fca9881f5901366..bb9683c39e90099aa445c8bd32949e1ce2010004 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C
@@ -163,11 +163,13 @@ JohnsonJacksonSchaeffer::nu
     const fvPatchList& patches = phase.mesh().boundary();
     const volVectorField& U = phase.U();
 
+    volScalarField::GeometricBoundaryField& nufBf = nuf.boundaryFieldRef();
+
     forAll(patches, patchi)
     {
         if (!patches[patchi].coupled())
         {
-            nuf.boundaryField()[patchi] =
+            nufBf[patchi] =
                 (
                     pf.boundaryField()[patchi]*sin(phi_.value())
                    /(
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index 8c072802d17df4132f48723f608d1ed8279bf223..0853ec845062bda1596700c50fbc24ca96dbfadc 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -152,11 +152,13 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
     const fvPatchList& patches = phase.mesh().boundary();
     const volVectorField& U = phase.U();
 
+    volScalarField::GeometricBoundaryField& nufBf = nuf.boundaryFieldRef();
+
     forAll(patches, patchi)
     {
         if (!patches[patchi].coupled())
         {
-            nuf.boundaryField()[patchi] =
+            nufBf[patchi] =
                 (
                     pf.boundaryField()[patchi]*sin(phi_.value())
                    /(
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index ff3be20ab10647ee6b6bc32793e374f06518b9e9..d0a5b023d1b23ac99f1c0a8ecdd9e207d6e8e118 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -295,7 +295,7 @@ Foam::RASModels::kineticTheoryModel::pPrime() const
     );
 
     volScalarField::GeometricBoundaryField& bpPrime =
-        tpPrime.ref().boundaryField();
+        tpPrime.ref().boundaryFieldRef();
 
     forAll(bpPrime, patchi)
     {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
index 9d3b53c782713b7657e257683c24b14570e7f09b..0edda3819d5a8dcded2a66d127a7a860e1c428b8 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
@@ -165,7 +165,7 @@ Foam::RASModels::phasePressureModel::pPrime() const
     );
 
     volScalarField::GeometricBoundaryField& bpPrime =
-        tpPrime.ref().boundaryField();
+        tpPrime.ref().boundaryFieldRef();
 
     forAll(bpPrime, patchi)
     {
@@ -193,7 +193,7 @@ Foam::RASModels::phasePressureModel::pPrimef() const
     );
 
    surfaceScalarField::GeometricBoundaryField& bpPrime =
-       tpPrime.ref().boundaryField();
+       tpPrime.ref().boundaryFieldRef();
 
     forAll(bpPrime, patchi)
     {
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 597d15e7ff04f702098af467fc8830cfde9eff43..fc86ede74b3155781281f6c9c98dfacd08c3112b 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -316,11 +316,13 @@ void Foam::twoPhaseSystem::solve()
             )
         );
 
+        surfaceScalarField::GeometricBoundaryField& alphaPhic1Bf =
+            alphaPhic1.boundaryFieldRef();
+
         // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhic1.boundaryField(), patchi)
+        forAll(alphaPhic1Bf, patchi)
         {
-            fvsPatchScalarField& alphaPhic1p =
-                alphaPhic1.boundaryField()[patchi];
+            fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi];
 
             if (!alphaPhic1p.coupled())
             {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
index b4188777437f2b2243d6ab8786c44adad822d8e4..eb3eeb562a55366e4e2980ab1356d4600b4eb0c5 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H
@@ -146,17 +146,25 @@ while (pimple.correct())
     surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99));
     surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar));
 
-    forAll(mesh.boundary(), patchi)
     {
-        // Set ddtPhiCorr to 0 on non-coupled boundaries
-        if
-        (
-            !mesh.boundary()[patchi].coupled()
-         || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
-        )
+        surfaceScalarField::GeometricBoundaryField& phiCorrCoeff1Bf =
+            phiCorrCoeff1.boundaryFieldRef();
+
+        surfaceScalarField::GeometricBoundaryField& phiCorrCoeff2Bf =
+            phiCorrCoeff2.boundaryFieldRef();
+
+        forAll(mesh.boundary(), patchi)
         {
-            phiCorrCoeff1.boundaryField()[patchi] = 0;
-            phiCorrCoeff2.boundaryField()[patchi] = 0;
+            // Set ddtPhiCorr to 0 on non-coupled boundaries
+            if
+            (
+               !mesh.boundary()[patchi].coupled()
+             || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
+            )
+            {
+                phiCorrCoeff1Bf[patchi] = 0;
+                phiCorrCoeff2Bf[patchi] = 0;
+            }
         }
     }
 
@@ -212,7 +220,7 @@ while (pimple.correct())
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
     (
-        p_rgh.boundaryField(),
+        p_rgh.boundaryFieldRef(),
         (
             phiHbyA.boundaryField()
           - (
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
index 6855b6fb962539e7c9a22c6565fe13f28a81e2e0..bb78b2799d43a3e750514d137ee02fa05433d203 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H
@@ -200,7 +200,7 @@ while (pimple.correct())
     // Update the fixedFluxPressure BCs to ensure flux consistency
     setSnGrad<fixedFluxPressureFvPatchScalarField>
     (
-        p_rgh.boundaryField(),
+        p_rgh.boundaryFieldRef(),
         (
             phiHbyA.boundaryField()
           - (
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index 6a5a474b3462459b792d79f07a004b0790417767..8d915bc39e9a9bb15a750e01d76ad4b2dd80876e 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -279,7 +279,7 @@ Foam::RASModels::kineticTheoryModel::pPrime() const
     );
 
     volScalarField::GeometricBoundaryField& bpPrime =
-        tpPrime.ref().boundaryField();
+        tpPrime.ref().boundaryFieldRef();
 
     forAll(bpPrime, patchi)
     {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
index 0c5906d905315af28178027aef65ed8a35bab062..3245dec4dbf7c82ceb9924f70e4adaffa8168836 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
@@ -171,7 +171,7 @@ Foam::RASModels::phasePressureModel::pPrime() const
     );
 
     volScalarField::GeometricBoundaryField& bpPrime =
-        tpPrime.ref().boundaryField();
+        tpPrime.ref().boundaryFieldRef();
 
     forAll(bpPrime, patchi)
     {
@@ -199,7 +199,7 @@ Foam::RASModels::phasePressureModel::pPrimef() const
     );
 
    surfaceScalarField::GeometricBoundaryField& bpPrime =
-       tpPrime.ref().boundaryField();
+       tpPrime.ref().boundaryFieldRef();
 
     forAll(bpPrime, patchi)
     {
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
index ecf674e6e0bfc04a049002207956a4942f0cea25..a3e187e7e37ae8193a53a2382d12780f73be6be8 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -36,6 +36,9 @@ void Foam::BlendedInterfacialModel<modelType>::correctFixedFluxBCs
     GeometricField& field
 ) const
 {
+    typename GeometricField::GeometricBoundaryField& fieldBf =
+        field.boundaryFieldRef();
+
     forAll(pair_.phase1().phi().boundaryField(), patchi)
     {
         if
@@ -46,7 +49,7 @@ void Foam::BlendedInterfacialModel<modelType>::correctFixedFluxBCs
             )
         )
         {
-            field.boundaryField()[patchi] = Zero;
+            fieldBf[patchi] = Zero;
         }
     }
 }
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 76af4a310a337ec58179b702556a7b86294f10cf..a8e810f73c88e6d58e72cb9b3d8c534b135f785c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -446,11 +446,13 @@ void Foam::twoPhaseSystem::solve()
             )
         );
 
+        surfaceScalarField::GeometricBoundaryField& alphaPhic1Bf =
+            alphaPhic1.boundaryFieldRef();
+
         // Ensure that the flux at inflow BCs is preserved
-        forAll(alphaPhic1.boundaryField(), patchi)
+        forAll(alphaPhic1Bf, patchi)
         {
-            fvsPatchScalarField& alphaPhic1p =
-                alphaPhic1.boundaryField()[patchi];
+            fvsPatchScalarField& alphaPhic1p = alphaPhic1Bf[patchi];
 
             if (!alphaPhic1p.coupled())
             {
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
index 2f1b71f7f1f4f52e6302deb28be271025f8a5854..1d2c3e825a17c83a99f7f33653499b2417a81992 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
@@ -429,30 +429,28 @@ public:
 
     // Member Functions
 
-        //- Return dimensioned internal field
+        //- Return a reference to the dimensioned internal field
+        //  Note: this increments the event counter and checks the
+        //  old-time fields; avoid in loops.
         DimensionedInternalField& dimensionedInternalField();
 
-        //- Return dimensioned internal field
+        //- Return a const-reference to the dimensioned internal field
         inline const DimensionedInternalField& dimensionedInternalField() const;
 
-        //- Return internal field
+        //- Return a reference to the internal field
+        //  Note: this increments the event counter and checks the
+        //  old-time fields; avoid in loops.
         InternalField& internalField();
 
-        //- Return internal field
+        //- Return a const-reference to the  internal field
         inline const InternalField& internalField() const;
 
-        //- Return reference to GeometricBoundaryField
+        //- Return a reference to the boundary field
+        //  Note: this increments the event counter and checks the
+        //  old-time fields; avoid in loops.
         GeometricBoundaryField& boundaryFieldRef();
 
-        //- Return reference to GeometricBoundaryField
-        #ifndef BOUNDARY_FIELD_REF
-        GeometricBoundaryField& boundaryField()
-        {
-            return boundaryFieldRef();
-        }
-        #endif
-
-        //- Return reference to GeometricBoundaryField for const field
+        //- Return const-reference to the boundary field
         inline const GeometricBoundaryField& boundaryField() const;
 
         //- Return the time index of the field
diff --git a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
index 32f44822474dd81f308d578f1356a67319c8da08..01b55287f15d36c4d441196f50c4de71bedfffa9 100644
--- a/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
+++ b/src/TurbulenceModels/phaseCompressible/RAS/mixtureKEpsilon/mixtureKEpsilon.C
@@ -196,7 +196,7 @@ void mixtureKEpsilon<BasicTurbulenceModel>::correctInletOutlet
     const volScalarField& refVsf
 ) const
 {
-    volScalarField::GeometricBoundaryField& bf = vsf.boundaryField();
+    volScalarField::GeometricBoundaryField& bf = vsf.boundaryFieldRef();
     const volScalarField::GeometricBoundaryField& refBf =
         refVsf.boundaryField();
 
diff --git a/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C b/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C
index 760a33546586571d47f35acaafa76690faa745ea..e9f1b40ecf0c41d0bfdbc4a7ccd41519b96f5015 100644
--- a/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C
+++ b/src/finiteVolume/cfdTools/general/constrainPressure/constrainPressure.C
@@ -44,7 +44,7 @@ void Foam::constrainPressure
 {
     const fvMesh& mesh = p.mesh();
 
-    volScalarField::GeometricBoundaryField& pBf = p.boundaryField();
+    volScalarField::GeometricBoundaryField& pBf = p.boundaryFieldRef();
 
     const volVectorField::GeometricBoundaryField& UBf = U.boundaryField();
     const surfaceScalarField::GeometricBoundaryField& phiHbyABf =
diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index 6673e5d983e2ea51fc69570874219f4d91cdbc39..b18852a1ae18862dedd12a1bcdae1ae00195780d 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -609,7 +609,7 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
         limitSum(phiPsiCorrsInternal);
     }
 
-    surfaceScalarField::GeometricBoundaryField& bfld =
+    const surfaceScalarField::GeometricBoundaryField& bfld =
         phiPsiCorrs[0].boundaryField();
 
     forAll(bfld, patchi)
@@ -622,7 +622,7 @@ void Foam::MULES::limitSum(SurfaceScalarFieldList& phiPsiCorrs)
                 phiPsiCorrsPatch.set
                 (
                     phasei,
-                    &phiPsiCorrs[phasei].boundaryField()[patchi]
+                    &phiPsiCorrs[phasei].boundaryFieldRef()[patchi]
                 );
             }