diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
index e003518593..c8927bfcf1 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C
@@ -391,12 +391,6 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
 }
 
 
-// Correction for the boundary condition on the unit normal nHat on
-// walls to produce the correct contact angle.
-
-// The dynamic contact angle is calculated from the component of the
-// velocity on the direction of the interface, parallel to the wall.
-
 void Foam::multiphaseMixture::correctContactAngle
 (
     const phase& alpha1,
@@ -404,21 +398,42 @@ void Foam::multiphaseMixture::correctContactAngle
     surfaceVectorField::Boundary& nHatb
 ) const
 {
-    const volScalarField::Boundary& gbf
-        = alpha1.boundaryField();
+    const volScalarField::Boundary& a1bf = alpha1.boundaryField();
+    const volScalarField::Boundary& a2bf = alpha2.boundaryField();
 
     const fvBoundaryMesh& boundary = mesh_.boundary();
 
     forAll(boundary, patchi)
     {
-        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
+        if
+        (
+            isA<alphaContactAngleFvPatchScalarField>(a1bf[patchi])
+         || isA<alphaContactAngleFvPatchScalarField>(a2bf[patchi])
+        )
         {
+            if
+            (
+                isA<alphaContactAngleFvPatchScalarField>(a1bf[patchi])
+             && isA<alphaContactAngleFvPatchScalarField>(a2bf[patchi])
+            )
+            {
+                FatalErrorInFunction
+                    << "alphaContactAngle boundary condition "
+                       "specified on patch " << boundary[patchi].name()
+                    << " for both " << alpha1.name() << " and " << alpha2.name()
+                    << nl << "which may be inconsistent."
+                    << exit(FatalError);
+            }
+
             const alphaContactAngleFvPatchScalarField& acap =
-                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
+                isA<alphaContactAngleFvPatchScalarField>(a1bf[patchi])
+              ? refCast<const alphaContactAngleFvPatchScalarField>(a1bf[patchi])
+              : refCast<const alphaContactAngleFvPatchScalarField>(a2bf[patchi])
+              ;
 
             vectorField& nHatPatch = nHatb[patchi];
 
-            vectorField AfHatPatch
+            const vectorField AfHatPatch
             (
                 mesh_.Sf().boundaryField()[patchi]
                /mesh_.magSf().boundaryField()[patchi]
@@ -437,18 +452,19 @@ void Foam::multiphaseMixture::correctContactAngle
                     << exit(FatalError);
             }
 
-            bool matched = (tp.key().first() == alpha1.name());
+            const bool matched = (tp.key().first() == alpha1.name());
+
+            const scalar theta0 = degToRad(tp().theta0(matched));
 
-            scalar theta0 = degToRad(tp().theta0(matched));
             scalarField theta(boundary[patchi].size(), theta0);
 
-            scalar uTheta = tp().uTheta();
+            const scalar uTheta = tp().uTheta();
 
             // Calculate the dynamic contact angle if required
             if (uTheta > small)
             {
-                scalar thetaA = degToRad(tp().thetaA(matched));
-                scalar thetaR = degToRad(tp().thetaR(matched));
+                const scalar thetaA = degToRad(tp().thetaA(matched));
+                const scalar thetaR = degToRad(tp().thetaR(matched));
 
                 // Calculated the component of the velocity parallel to the wall
                 vectorField Uwall
@@ -469,7 +485,7 @@ void Foam::multiphaseMixture::correctContactAngle
 
                 // Calculate Uwall resolved normal to the interface parallel to
                 // the interface
-                scalarField uwall(nWall & Uwall);
+                const scalarField uwall(nWall & Uwall);
 
                 theta += (thetaA - thetaR)*tanh(uwall/uTheta);
             }
@@ -477,9 +493,9 @@ void Foam::multiphaseMixture::correctContactAngle
 
             // Reset nHatPatch to correspond to the contact angle
 
-            scalarField a12(nHatPatch & AfHatPatch);
+            const scalarField a12(nHatPatch & AfHatPatch);
 
-            scalarField b1(cos(theta));
+            const scalarField b1(cos(theta));
 
             scalarField b2(nHatPatch.size());
 
@@ -488,10 +504,10 @@ void Foam::multiphaseMixture::correctContactAngle
                 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
             }
 
-            scalarField det(1.0 - a12*a12);
+            const scalarField det(1.0 - a12*a12);
 
-            scalarField a((b1 - a12*b2)/det);
-            scalarField b((b2 - a12*b1)/det);
+            const scalarField a((b1 - a12*b2)/det);
+            const scalarField b((b2 - a12*b1)/det);
 
             nHatPatch = a*AfHatPatch + b*nHatPatch;
 
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
index b28040a8c9..e85f7c38a1 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.H
@@ -173,6 +173,11 @@ private:
             const volScalarField& alpha2
         ) const;
 
+        //- Correction for the boundary condition on the unit normal nHat on
+        // walls to produce the correct contact angle.
+        //
+        // The dynamic contact angle is calculated from the component of the
+        // velocity on the direction of the interface, parallel to the wall.
         void correctContactAngle
         (
             const phase& alpha1,