diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/contErrs.H b/applications/solvers/multiphase/twoPhaseEulerFoam/contErrs.H
new file mode 100644
index 0000000000000000000000000000000000000000..83ca7cea924a191d36cba7728794ac8f38a34cb4
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/contErrs.H
@@ -0,0 +1,13 @@
+volScalarField contErr1
+(
+    "contErr1",
+    fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+  - (fvOptions(alpha1, rho1)&rho1)
+);
+
+volScalarField contErr2
+(
+    "contErr2",
+    fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+  - (fvOptions(alpha2, rho2)&rho2)
+);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/correctContErrs.H b/applications/solvers/multiphase/twoPhaseEulerFoam/correctContErrs.H
new file mode 100644
index 0000000000000000000000000000000000000000..75355b3963fb4cca83dcc5fcdc1d843043a88a46
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/correctContErrs.H
@@ -0,0 +1,7 @@
+contErr1 =
+    fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+  - (fvOptions(alpha1, rho1)&rho1);
+
+contErr2 =
+    fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+  - (fvOptions(alpha2, rho2)&rho2);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index edac7acd49f544b2895a4d2c045617cc250d5012..0924c5d178c86852322b55870127a939bd1e1da8 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -10,6 +10,9 @@ surfaceScalarField alpharAU2f(fvc::interpolate(alpha2*rAU2));
 // --- Pressure corrector loop
 while (pimple.correct())
 {
+    // Update continuity errors due to temperature changes
+    #include "correctContErrs.H"
+
     volVectorField HbyA1
     (
         IOobject::groupName("HbyA", phase1.name()),
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index dd3a8fb710fd25458663cfd556d7275c9ff7d651..f9effa437ba8482ba82094b3838a75b0544dc1f0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -76,20 +76,7 @@ int main(int argc, char *argv[])
             fluid.solve();
             fluid.correct();
 
-            volScalarField contErr1
-            (
-                "contErr1",
-                fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
-              - (fvOptions(alpha1, rho1)&rho1)
-            );
-
-            volScalarField contErr2
-            (
-                "contErr2",
-                fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
-              - (fvOptions(alpha2, rho2)&rho2)
-            );
-
+            #include "contErrs.H"
             #include "UEqns.H"
             #include "EEqns.H"
             #include "pEqn.H"
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/T.particles b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/T.particles
index 58313dfb2b5c6f5d87e45ee056240033c9e61406..5762b51bd8e92a7a48f6b3f63a39d7192a0b8e0c 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/T.particles
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/0/T.particles
@@ -27,7 +27,10 @@ boundaryField
 
     outlet
     {
-        type               zeroGradient;
+        type               inletOutlet;
+        phi                phi.particles;
+        inletValue         uniform 300;
+        value              uniform 300;
     }
 
     walls
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/T.particles b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/T.particles
index 58313dfb2b5c6f5d87e45ee056240033c9e61406..5762b51bd8e92a7a48f6b3f63a39d7192a0b8e0c 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/T.particles
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/0/T.particles
@@ -27,7 +27,10 @@ boundaryField
 
     outlet
     {
-        type               zeroGradient;
+        type               inletOutlet;
+        phi                phi.particles;
+        inletValue         uniform 300;
+        value              uniform 300;
     }
 
     walls