diff --git a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C
index 02290c7038..348c3de101 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C
@@ -43,7 +43,7 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
 :
     rhoThermo::composite(U.mesh(), word::null),
     twoPhaseMixture(U.mesh(), *this),
-    interfaceProperties(alpha1(), U, *this),
+    interfaceProperties(alpha1(), alpha2(), U, *this),
     thermo1_(nullptr),
     thermo2_(nullptr),
     Alpha1_
diff --git a/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C b/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C
index d009acbef9..863c2a80f6 100644
--- a/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C
+++ b/src/twoPhaseModels/immiscibleIncompressibleTwoPhaseMixture/immiscibleIncompressibleTwoPhaseMixture.C
@@ -36,7 +36,7 @@ immiscibleIncompressibleTwoPhaseMixture
 )
 :
     incompressibleTwoPhaseMixture(U, phi),
-    interfaceProperties(alpha1(), U, *this)
+    interfaceProperties(alpha1(), alpha2(), U, *this)
 {}
 
 
diff --git a/src/twoPhaseModels/interfaceProperties/interfaceProperties.C b/src/twoPhaseModels/interfaceProperties/interfaceProperties.C
index 8890e71ff4..4a1175ce01 100644
--- a/src/twoPhaseModels/interfaceProperties/interfaceProperties.C
+++ b/src/twoPhaseModels/interfaceProperties/interfaceProperties.C
@@ -43,30 +43,28 @@ void Foam::interfaceProperties::correctContactAngle
 (
     surfaceVectorField::Boundary& nHatb,
     const surfaceVectorField::Boundary& gradAlphaf
-) const
+)
 {
     const fvMesh& mesh = alpha1_.mesh();
-    const volScalarField::Boundary& abf = alpha1_.boundaryField();
+    volScalarField::Boundary& a1bf = alpha1_.boundaryFieldRef();
+    volScalarField::Boundary& a2bf = alpha2_.boundaryFieldRef();
 
     const fvBoundaryMesh& boundary = mesh.boundary();
 
     forAll(boundary, patchi)
     {
-        if (isA<alphaContactAngleFvPatchScalarField>(abf[patchi]))
+        if (isA<alphaContactAngleFvPatchScalarField>(a1bf[patchi]))
         {
-            alphaContactAngleFvPatchScalarField& acap =
-                const_cast<alphaContactAngleFvPatchScalarField&>
+            alphaContactAngleFvPatchScalarField& a1cap =
+                refCast<alphaContactAngleFvPatchScalarField>
                 (
-                    refCast<const alphaContactAngleFvPatchScalarField>
-                    (
-                        abf[patchi]
-                    )
+                    a1bf[patchi]
                 );
 
             fvsPatchVectorField& nHatp = nHatb[patchi];
             const scalarField theta
             (
-                degToRad(acap.theta(U_.boundaryField()[patchi], nHatp))
+                degToRad(a1cap.theta(U_.boundaryField()[patchi], nHatp))
             );
 
             const vectorField nf
@@ -93,8 +91,9 @@ void Foam::interfaceProperties::correctContactAngle
             nHatp = a*nf + b*nHatp;
             nHatp /= (mag(nHatp) + deltaN_.value());
 
-            acap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
-            acap.evaluate();
+            a1cap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
+            a1cap.evaluate();
+            a2bf[patchi] = 1 - a1cap;
         }
     }
 }
@@ -148,7 +147,8 @@ void Foam::interfaceProperties::calculateK()
 
 Foam::interfaceProperties::interfaceProperties
 (
-    const volScalarField& alpha1,
+    volScalarField& alpha1,
+    volScalarField& alpha2,
     const volVectorField& U,
     const IOdictionary& dict
 )
@@ -164,6 +164,7 @@ Foam::interfaceProperties::interfaceProperties
     ),
 
     alpha1_(alpha1),
+    alpha2_(alpha2),
     U_(U),
 
     nHatf_
diff --git a/src/twoPhaseModels/interfaceProperties/interfaceProperties.H b/src/twoPhaseModels/interfaceProperties/interfaceProperties.H
index 1f06363f80..429a058228 100644
--- a/src/twoPhaseModels/interfaceProperties/interfaceProperties.H
+++ b/src/twoPhaseModels/interfaceProperties/interfaceProperties.H
@@ -66,7 +66,8 @@ class interfaceProperties
         //- Stabilisation for normalisation of the interface normal
         const dimensionedScalar deltaN_;
 
-        const volScalarField& alpha1_;
+        volScalarField& alpha1_;
+        volScalarField& alpha2_;
         const volVectorField& U_;
         surfaceScalarField nHatf_;
         volScalarField K_;
@@ -81,7 +82,7 @@ class interfaceProperties
         (
             surfaceVectorField::Boundary& nHat,
             const surfaceVectorField::Boundary& gradAlphaf
-        ) const;
+        );
 
         //- Re-calculate the interface curvature
         void calculateK();
@@ -94,7 +95,8 @@ public:
         //- Construct from volume fraction field gamma and IOdictionary
         interfaceProperties
         (
-            const volScalarField& alpha1,
+            volScalarField& alpha1,
+            volScalarField& alpha2,
             const volVectorField& U,
             const IOdictionary&
         );