diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index a4e9e0b40f6c77500d7b55adad6daa8d399c3279..3c4cbb0d32e48f0ff38415408ec669119b4f9f03 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -4,8 +4,8 @@
         p =
         (
             rho
-          - (1.0 - gamma)*rhol0
-          - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+          - gamma2*rhol0
+          - ((gamma*psiv + gamma2*psil) - psi)*pSat
         )/psi;
     }
 
@@ -57,8 +57,8 @@
     p =
     (
         rho
-      - (1.0 - gamma)*rhol0
-      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+      - gamma2*rhol0
+      - ((gamma*psiv + gamma2*psil) - psi)*pSat
     )/psi;
 
     p.correctBoundaryConditions();
diff --git a/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H b/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H
index c5e93c7efb3d537892e22228143b8103e9ce3a21..ce618ee40f2c3a04bf660d08c385389ecdebd570 100644
--- a/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H
+++ b/applications/solvers/multiphase/cavitatingFoam/continuityErrs.H
@@ -1,5 +1,5 @@
 {
-    volScalarField thermoRho = psi*p + (1.0 - gamma)*rhol0;
+    volScalarField thermoRho = psi*p + gamma2*rhol0;
 
     dimensionedScalar totalMass = fvc::domainIntegrate(rho);
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/createFields.H b/applications/solvers/multiphase/cavitatingFoam/createFields.H
index dbacf1dbd3f1c416bffccc6e2c3b5e7f67d4baf0..c7241406504a6d32d08a9c28b552d0fcdcc3eb44 100644
--- a/applications/solvers/multiphase/cavitatingFoam/createFields.H
+++ b/applications/solvers/multiphase/cavitatingFoam/createFields.H
@@ -49,6 +49,8 @@
     volScalarField& gamma(twoPhaseProperties.alpha1());
     gamma.oldTime();
 
+    volScalarField& gamma2(twoPhaseProperties.alpha2());
+
     Info<< "Creating compressibilityModel\n" << endl;
     autoPtr<barotropicCompressibilityModel> psiModel =
         barotropicCompressibilityModel::New
@@ -62,8 +64,8 @@
     rho == max
     (
         psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+      + gamma2*rhol0
+      + ((gamma*psiv + gamma2*psil) - psi)*pSat,
         rhoMin
     );
 
diff --git a/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H b/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H
index b259ddd3222ff00fe49a83854a1bcae1e84ba240..4edda7336de09a5e212dbded85bcc70c2ddee6dc 100644
--- a/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H
+++ b/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H
@@ -1,5 +1,6 @@
 {
     gamma = max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));
+    gamma2 = 1.0 - gamma;
 
     Info<< "max-min gamma: " << max(gamma).value()
         << " " << min(gamma).value() << endl;
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index ea35f79fa14c505d87f62be95f2bab9fd1c56bf9..b2c7d953b5368b2b46adb91223c06d9ceab4242d 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -4,8 +4,8 @@
         p =
         (
             rho
-          - (1.0 - gamma)*rhol0
-          - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+          - gamma2*rhol0
+          - ((gamma*psiv + gamma2*psil) - psi)*pSat
         )/psi;
     }
 
@@ -49,8 +49,8 @@
     rho == max
     (
         psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+      + gamma2*rhol0
+      + ((gamma*psiv + gamma2*psil) - psi)*pSat,
         rhoMin
     );
 
@@ -59,8 +59,8 @@
     p =
     (
         rho
-      - (1.0 - gamma)*rhol0
-      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+      - gamma2*rhol0
+      - ((gamma*psiv + gamma2*psil) - psi)*pSat
     )/psi;
 
     p.correctBoundaryConditions();
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
index 5f5ac824b8aacfdd77227f03bdcba07026870939..faae19767035025223dbd22b6be483d7dc6d6ace 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
@@ -18,7 +18,7 @@
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha1,
                 alpharScheme
             )
@@ -27,6 +27,7 @@
         MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
         //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
 
+        alpha2 = 1.0 - alpha1;
         rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
     }
 
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
index 6c82f940499ea39567db1e400274832dfa354b0f..57c78027a471f712cf06e17486b62fcd05cd3014 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@@ -23,4 +23,4 @@ else
     #include "alphaEqn.H"
 }
 
-rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho == alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index b96dcf898daaeaec279adcce7ac3f5916623089a..a2720e20eeec598b64ac4e783e175670069f98da 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -18,7 +18,7 @@
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha1,
                 alpharScheme
             )
@@ -26,6 +26,7 @@
 
         MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
 
+        alpha2 = 1.0 - alpha1;
         rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
     }
 
diff --git a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
index 6c82f940499ea39567db1e400274832dfa354b0f..57c78027a471f712cf06e17486b62fcd05cd3014 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H
@@ -23,4 +23,4 @@ else
     #include "alphaEqn.H"
 }
 
-rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho == alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/interFoam/createFields.H b/applications/solvers/multiphase/interFoam/createFields.H
index df1edf04e1e4c1ef7f1cf802d4645bbe2c0ca122..c607c6854891c7b83498075d75ee734abe0cb118 100644
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@@ -33,6 +33,7 @@
     twoPhaseMixture twoPhaseProperties(U, phi);
 
     volScalarField& alpha1(twoPhaseProperties.alpha1());
+    volScalarField& alpha2(twoPhaseProperties.alpha2());
 
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
@@ -48,7 +49,7 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
+        alpha1*rho1 + alpha2*rho2,
         alpha1.boundaryField().types()
     );
     rho.oldTime();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
index 139646e43c6023c1486e3ee53dfdb0a14c63e158..6bc5953a4a6424f3bd28406e6f09ed88a8a30680 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaDiffusionEqn.H
@@ -12,7 +12,8 @@
 
     alpha1Eqn.solve();
 
+    alpha2 = 1.0 - alpha1;
     rhoPhi += alpha1Eqn.flux()*(rho1 - rho2);
 }
 
-rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho = alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
index 6c82f940499ea39567db1e400274832dfa354b0f..57c78027a471f712cf06e17486b62fcd05cd3014 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H
@@ -23,4 +23,4 @@ else
     #include "alphaEqn.H"
 }
 
-rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+rho == alpha1*rho1 + alpha2*rho2;
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
index 0d01b9a9e5ed82abc4a0420ece0ca464830b7b97..e98f06e74e590fb3a2b8a635fc28324130967251 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@@ -32,6 +32,7 @@
     twoPhaseMixture twoPhaseProperties(U, phi);
 
     volScalarField& alpha1(twoPhaseProperties.alpha1());
+    volScalarField& alpha2(twoPhaseProperties.alpha2());
 
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
@@ -42,7 +43,7 @@
     dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab"));
 
     // Need to store rho for ddt(rho, U)
-    volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2);
+    volScalarField rho("rho", alpha1*rho1 + alpha2*rho2);
     rho.oldTime();
 
 
diff --git a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C
index 4ec3ce72f85dc8b19a6392c4111eaf05b0247dfb..4b698138a77cb06e70294325373311ab7c19cee3 100644
--- a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C
+++ b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.C
@@ -54,7 +54,8 @@ Foam::twoPhaseMixture::twoPhaseMixture
 (
     const volVectorField& U,
     const surfaceScalarField& phi,
-    const word& alpha1Name
+    const word& alpha1Name,
+    const word& alpha2Name
 )
 :
     transportModel(U, phi),
@@ -102,6 +103,17 @@ Foam::twoPhaseMixture::twoPhaseMixture
         U_.mesh()
     ),
 
+    alpha2_
+    (
+        IOobject
+        (
+            found("phases") ? word("alpha" + phase2Name_) : alpha2Name,
+            U_.time().timeName(),
+            U_.db()
+        ),
+        1.0 - alpha1_
+    ),
+
     nu_
     (
         IOobject
diff --git a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H
index 025d95f7011b609c3fd3bb0bc914216b74bdaeb4..265f862bce2c39345182e1087c625b211ba24b52 100644
--- a/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H
+++ b/src/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H
@@ -70,6 +70,7 @@ protected:
         const surfaceScalarField& phi_;
 
         volScalarField alpha1_;
+        volScalarField alpha2_;
 
         volScalarField nu_;
 
@@ -89,7 +90,8 @@ public:
         (
             const volVectorField& U,
             const surfaceScalarField& phi,
-            const word& alpha1Name = "alpha1"
+            const word& alpha1Name = "alpha1",
+            const word& alpha2Name = "alpha2"
         );
 
 
@@ -122,6 +124,18 @@ public:
             return alpha1_;
         }
 
+        //- Return the phase-fraction of phase 2
+        const volScalarField& alpha2() const
+        {
+            return alpha2_;
+        }
+
+        //- Return the phase-fraction of phase 2
+        volScalarField& alpha2()
+        {
+            return alpha2_;
+        }
+
         //- Return const-access to phase1 viscosityModel
         const viscosityModel& nuModel1() const
         {