diff --git a/applications/solvers/modules/Allwmake b/applications/solvers/modules/Allwmake
index 9b2cb328a7..7d0d8dc874 100755
--- a/applications/solvers/modules/Allwmake
+++ b/applications/solvers/modules/Allwmake
@@ -16,6 +16,7 @@ incompressibleVoF/Allwmake $targetType $*
 compressibleVoF/Allwmake $targetType $*
 wmake $targetType multiphaseVoFSolver
 wmake $targetType incompressibleMultiphaseVoF
+wmake $targetType compressibleMultiphaseVoF
 multiphaseEuler/Allwmake $targetType $*
 wmake $targetType solid
 solidDisplacement/Allwmake $targetType $*
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/Make/files b/applications/solvers/modules/compressibleMultiphaseVoF/Make/files
new file mode 100644
index 0000000000..81da7fe42a
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/Make/files
@@ -0,0 +1,11 @@
+compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.C
+compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixtureThermo/compressibleMultiphaseVoFMixtureThermo.C
+compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C
+
+alphaPredictor.C
+momentumPredictor.C
+thermophysicalPredictor.C
+pressureCorrector.C
+compressibleMultiphaseVoF.C
+
+LIB = $(FOAM_LIBBIN)/libcompressibleMultiphaseVoF
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/Make/options b/applications/solvers/modules/compressibleMultiphaseVoF/Make/options
new file mode 100644
index 0000000000..8073e24c1a
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/Make/options
@@ -0,0 +1,28 @@
+EXE_INC = \
+    -I$(FOAM_SOLVERS)/modules/multiphaseVoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/VoFSolver/lnInclude \
+    -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
+    -I$(LIB_SRC)/physicalProperties/lnInclude \
+    -I$(LIB_SRC)/multiphaseModels/multiphaseProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/phaseIncompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lmultiphaseVoFSolver \
+    -lphysicalProperties \
+    -lmultiphaseProperties \
+    -lfluidThermophysicalModels \
+    -linterfaceCompression \
+    -lmomentumTransportModels \
+    -lcompressibleMomentumTransportModels \
+    -lphaseIncompressibleMomentumTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvModels \
+    -lfvConstraints \
+    -lsampling
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/alphaPredictor.C b/applications/solvers/modules/compressibleMultiphaseVoF/alphaPredictor.C
new file mode 100644
index 0000000000..b388a8bdd4
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/alphaPredictor.C
@@ -0,0 +1,277 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleMultiphaseVoF.H"
+#include "subCycle.H"
+#include "CMULES.H"
+#include "fvcFlux.H"
+#include "fvcMeshPhi.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleMultiphaseVoF::alphaSolve
+(
+    const dictionary& alphaControls
+)
+{
+    const scalar cAlpha(alphaControls.lookup<scalar>("cAlpha"));
+
+    const word alphaScheme("div(phi,alpha)");
+    const word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phic(mag(phi/mesh.magSf()));
+    phic = min(cAlpha*phic, max(phic));
+
+    UPtrList<const volScalarField> alphas(phases.size());
+    PtrList<surfaceScalarField> alphaPhis(phases.size());
+
+    forAll(phases, phasei)
+    {
+        const compressibleVoFphase& alpha = phases[phasei];
+
+        alphas.set(phasei, &alpha);
+
+        alphaPhis.set
+        (
+            phasei,
+            new surfaceScalarField
+            (
+                "phi" + alpha.name() + "Corr",
+                fvc::flux
+                (
+                    phi,
+                    alpha,
+                    alphaScheme
+                )
+            )
+        );
+
+        surfaceScalarField& alphaPhi = alphaPhis[phasei];
+
+        forAll(phases, phasej)
+        {
+            compressibleVoFphase& alpha2 = phases[phasej];
+
+            if (&alpha2 == &alpha) continue;
+
+            surfaceScalarField phir(phic*mixture.nHatf(alpha, alpha2));
+
+            alphaPhi += fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha,
+                alpharScheme
+            );
+        }
+
+        // Limit alphaPhi for each phase
+        MULES::limit
+        (
+            1.0/mesh.time().deltaT().value(),
+            geometricOneField(),
+            alpha,
+            phi,
+            alphaPhi,
+            zeroField(),
+            zeroField(),
+            oneField(),
+            zeroField(),
+            false
+        );
+    }
+
+    MULES::limitSum(alphas, alphaPhis, phi);
+
+    rhoPhi = Zero;
+
+    volScalarField sumAlpha
+    (
+        IOobject
+        (
+            "sumAlpha",
+            mesh.time().name(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar(dimless, 0)
+    );
+
+    const volScalarField divU(fvc::div(fvc::absolute(phi, U)));
+
+    forAll(phases, phasei)
+    {
+        compressibleVoFphase& alpha = phases[phasei];
+
+        surfaceScalarField& alphaPhi = alphaPhis[phasei];
+
+        volScalarField::Internal Sp
+        (
+            IOobject
+            (
+                "Sp",
+                mesh.time().name(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(alpha.dgdt().dimensions(), 0)
+        );
+
+        volScalarField::Internal Su
+        (
+            IOobject
+            (
+                "Su",
+                mesh.time().name(),
+                mesh
+            ),
+            // Divergence term is handled explicitly to be
+            // consistent with the explicit transport solution
+            divU.v()*min(alpha.v(), scalar(1))
+        );
+
+        {
+            const scalarField& dgdt = alpha.dgdt();
+
+            forAll(dgdt, celli)
+            {
+                if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
+                {
+                    Sp[celli] += dgdt[celli]*alpha[celli];
+                    Su[celli] -= dgdt[celli]*alpha[celli];
+                }
+                else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
+                {
+                    Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
+                }
+            }
+        }
+
+
+        forAll(phases, phasej)
+        {
+            const compressibleVoFphase& alpha2 = phases[phasej];
+
+            if (&alpha2 == &alpha) continue;
+
+            const scalarField& dgdt2 = alpha2.dgdt();
+
+            forAll(dgdt2, celli)
+            {
+                if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
+                {
+                    Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
+                    Su[celli] += dgdt2[celli]*alpha[celli];
+                }
+                else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
+                {
+                    Sp[celli] += dgdt2[celli]*alpha2[celli];
+                }
+            }
+        }
+
+        MULES::explicitSolve
+        (
+            geometricOneField(),
+            alpha,
+            alphaPhi,
+            Sp,
+            Su
+        );
+
+        rhoPhi += fvc::interpolate(alpha.thermo().rho())*alphaPhi;
+
+        Info<< alpha.name() << " volume fraction, min, max = "
+            << alpha.weightedAverage(mesh.V()).value()
+            << ' ' << min(alpha).value()
+            << ' ' << max(alpha).value()
+            << endl;
+
+        sumAlpha += alpha;
+    }
+
+    // Correct the sum of the phase-fractions to avoid 'drift'
+    const volScalarField sumCorr(1.0 - sumAlpha);
+    forAll(phases, phasei)
+    {
+        compressibleVoFphase& alpha = phases[phasei];
+        alpha += alpha*sumCorr;
+    }
+}
+
+
+void Foam::solvers::compressibleMultiphaseVoF::alphaPredictor()
+{
+    const dictionary& alphaControls = mesh.solution().solverDict("alpha");
+
+    const label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
+
+    if (nAlphaSubCycles > 1)
+    {
+        surfaceScalarField rhoPhiSum
+        (
+            IOobject
+            (
+                "rhoPhiSum",
+                runTime.name(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar(rhoPhi.dimensions(), 0)
+        );
+
+        const dimensionedScalar totalDeltaT = runTime.deltaT();
+
+        List<volScalarField*> alphaPtrs(phases.size());
+        forAll(phases, phasei)
+        {
+            alphaPtrs[phasei] = &phases[phasei];
+        }
+
+        for
+        (
+            subCycle<volScalarField, subCycleFields> alphaSubCycle
+            (
+                alphaPtrs,
+                nAlphaSubCycles
+            );
+            !(++alphaSubCycle).end();
+        )
+        {
+            alphaSolve(alphaControls);
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+        }
+
+        rhoPhi = rhoPhiSum;
+    }
+    else
+    {
+        alphaSolve(alphaControls);
+    }
+
+    mixture.correct();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C
new file mode 100644
index 0000000000..5b00d558b7
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleMultiphaseVoF.H"
+#include "CorrectPhi.H"
+#include "geometricZeroField.H"
+#include "fvcDdt.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+    defineTypeNameAndDebug(compressibleMultiphaseVoF, 0);
+    addToRunTimeSelectionTable(solver, compressibleMultiphaseVoF, fvMesh);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solvers::compressibleMultiphaseVoF::compressibleMultiphaseVoF
+(
+    fvMesh& mesh
+)
+:
+    multiphaseVoFSolver
+    (
+        mesh,
+        autoPtr<multiphaseVoFMixture>
+        (
+            new compressibleMultiphaseVoFMixture(mesh)
+        )
+    ),
+
+    mixture
+    (
+        refCast<compressibleMultiphaseVoFMixture>
+        (
+            multiphaseVoFSolver::mixture
+        )
+    ),
+
+    phases(mixture.phases()),
+
+    p(mixture.p()),
+
+    pressureReference_
+    (
+        p,
+        p_rgh,
+        pimple.dict(),
+        false
+    ),
+
+    pMin
+    (
+        "pMin",
+        dimPressure,
+        mixture
+    ),
+
+    K("K", 0.5*magSqr(U)),
+
+    momentumTransport_
+    (
+        compressible::momentumTransportModel::New
+        (
+            rho,
+            U,
+            rhoPhi,
+            mixture
+        )
+    ),
+
+    momentumTransport(momentumTransport_())
+{
+    // Read the controls
+    read();
+
+    if (correctPhi)
+    {
+        rAU = new volScalarField
+        (
+            IOobject
+            (
+                "rAU",
+                runTime.name(),
+                mesh,
+                IOobject::READ_IF_PRESENT,
+                IOobject::AUTO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimTime/dimDensity, 1)
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solvers::compressibleMultiphaseVoF::~compressibleMultiphaseVoF()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleMultiphaseVoF::prePredictor()
+{
+    multiphaseVoFSolver::prePredictor();
+
+    contErr = fvc::ddt(rho)()() + fvc::div(rhoPhi)()();
+
+    forAll(mixture.phases(), phasei)
+    {
+        const volScalarField& rho = phases[phasei].thermo().rho();
+        contErr.ref() -= fvModels().source(phases[phasei], rho)&rho;
+    }
+
+    if (pimple.predictTransport())
+    {
+        momentumTransport.predict();
+    }
+}
+
+
+void Foam::solvers::compressibleMultiphaseVoF::postCorrector()
+{
+    if (pimple.correctTransport())
+    {
+        momentumTransport.correct();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H
new file mode 100644
index 0000000000..bbb0d54cb2
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H
@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::solvers::compressibleMultiphaseVoF
+
+Description
+    Solver module for the solution of multiple compressible, isothermal
+    immiscible fluids using a VOF (volume of fluid) phase-fraction based
+    interface capturing approach, with optional mesh motion and mesh topology
+    changes including adaptive re-meshing.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    A mixture approach for momentum transport is provided in which a single
+    laminar, RAS or LES model is selected to model the momentum stress.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient and steady simulations.
+
+SourceFiles
+    compressibleMultiphaseVoF.C
+
+See also
+    Foam::solvers::VoFSolver
+    Foam::solvers::multiphaseVoFSolver
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleMultiphaseVoF_H
+#define compressibleMultiphaseVoF_H
+
+#include "multiphaseVoFSolver.H"
+#include "compressibleMultiphaseVoFMixture.H"
+#include "compressibleMomentumTransportModels.H"
+#include "pressureReference.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solvers
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class compressibleMultiphaseVoF Declaration
+\*---------------------------------------------------------------------------*/
+
+class compressibleMultiphaseVoF
+:
+    public multiphaseVoFSolver
+{
+
+protected:
+
+    // Phase properties
+
+        //- The compressible two-phase mixture
+        compressibleMultiphaseVoFMixture& mixture;
+
+        //- Reference to the phases
+        UPtrListDictionary<compressibleVoFphase>& phases;
+
+
+    // Thermophysical properties
+
+        //- Reference to the mixture static pressure field
+        volScalarField& p;
+
+
+    // Pressure reference
+
+        //- Pressure reference
+        Foam::pressureReference pressureReference_;
+
+        //- Minimum pressure
+        dimensionedScalar pMin;
+
+
+    // Kinematic properties
+
+        //- Kinetic energy field
+        //  Used in the energy equation
+        volScalarField K;
+
+
+    // Momentum transport
+
+        //- Momentum transport model pointer
+        autoPtr<compressible::momentumTransportModel> momentumTransport_;
+
+        //- Momentum transport model
+        compressible::momentumTransportModel& momentumTransport;
+
+
+    // Cached temporary fields
+
+        //- Continuity error
+        tmp<volScalarField::Internal> contErr;
+
+
+    // Protected Member Functions
+
+        //- Return the pressure reference
+        virtual const Foam::pressureReference& pressureReference() const
+        {
+            return pressureReference_;
+        }
+
+        //- Compressible flow is divergent
+        virtual bool divergent()
+        {
+            return true;
+        }
+
+        //- Return the momentum equation stress term
+        virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U);
+
+        //- Solve for the phase-fractions
+        virtual void alphaPredictor();
+
+
+private:
+
+    // Private Member Functions
+
+        //- Solve for the phase-fractions
+        void alphaSolve(const dictionary& alphaControls);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("compressibleMultiphaseVoF");
+
+
+    // Constructors
+
+        //- Construct from region mesh
+        compressibleMultiphaseVoF(fvMesh& mesh);
+
+        //- Disallow default bitwise copy construction
+        compressibleMultiphaseVoF
+        (
+            const compressibleMultiphaseVoF&
+        ) = delete;
+
+
+    //- Destructor
+    virtual ~compressibleMultiphaseVoF();
+
+
+    // Member Functions
+
+        //- Called at the start of the PIMPLE loop
+        virtual void prePredictor();
+
+        //- Construct and optionally solve the momentum equation
+        virtual void momentumPredictor();
+
+        //- Construct and solve the energy equation,
+        //  convert to temperature
+        //  and update thermophysical and transport properties
+        virtual void thermophysicalPredictor();
+
+        //- Construct and solve the pressure equation in the PISO loop
+        virtual void pressureCorrector();
+
+        //- Correct the momentum and thermophysical transport modelling
+        virtual void postCorrector();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const compressibleMultiphaseVoF&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C
new file mode 100644
index 0000000000..be4046e187
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C
@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleMultiphaseVoFMixture.H"
+#include "surfaceInterpolate.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(compressibleMultiphaseVoFMixture, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::compressibleMultiphaseVoFMixture::compressibleMultiphaseVoFMixture
+(
+    const fvMesh& mesh
+)
+:
+    compressibleMultiphaseVoFMixtureThermo(mesh),
+
+    multiphaseVoFMixture(mesh, compressibleVoFphase::iNew(mesh, T())),
+
+    phases_(multiphaseVoFMixture::phases().convert<compressibleVoFphase>()),
+
+    rho_
+    (
+        IOobject
+        (
+            "rho",
+            mesh.time().name(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimDensity, 0)
+    )
+{
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::compressibleMultiphaseVoFMixture::nu() const
+{
+    volScalarField mu(phases_[0].Alpha()*phases_[0].thermo().mu());
+
+    for (label phasei=1; phasei<phases_.size(); phasei++)
+    {
+        mu += phases_[phasei].Alpha()*phases_[phasei].thermo().mu();
+    }
+
+    return mu/rho_;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::compressibleMultiphaseVoFMixture::nu
+(
+    const label patchi
+) const
+{
+    scalarField mu
+    (
+        phases_[0].Alpha().boundaryField()[patchi]
+       *phases_[0].thermo().mu(patchi)
+    );
+
+    for (label phasei=1; phasei<phases_.size(); phasei++)
+    {
+        mu +=
+            phases_[phasei].Alpha().boundaryField()[patchi]
+           *phases_[phasei].thermo().mu(patchi);
+    }
+
+    return mu/rho_.boundaryField()[patchi];
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::compressibleMultiphaseVoFMixture::alphaEff
+(
+    const volScalarField& nut
+) const
+{
+    tmp<volScalarField> talphaEff
+    (
+        phases_[0]
+       *(
+           phases_[0].thermo().kappa()
+         + phases_[0].thermo().rho()*phases_[0].thermo().Cp()*nut
+        )/phases_[0].thermo().Cv()
+    );
+
+    for (label phasei=1; phasei<phases_.size(); phasei++)
+    {
+        talphaEff.ref() +=
+            phases_[phasei]
+           *(
+               phases_[phasei].thermo().kappa()
+             + phases_[phasei].thermo().rho()*phases_[phasei].thermo().Cp()*nut
+            )/phases_[phasei].thermo().Cv();
+    }
+
+    return talphaEff;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::compressibleMultiphaseVoFMixture::rCv() const
+{
+    tmp<volScalarField> trCv(phases_[0]/phases_[0].thermo().Cv());
+
+    for (label phasei=1; phasei<phases_.size(); phasei++)
+    {
+        trCv.ref() += phases_[phasei]/phases_[phasei].thermo().Cv();
+    }
+
+    return trCv;
+}
+
+
+void Foam::compressibleMultiphaseVoFMixture::correctThermo()
+{
+    forAll(phases_, phasei)
+    {
+        phases_[phasei].correct(p(), T());
+    }
+}
+
+
+void Foam::compressibleMultiphaseVoFMixture::correct()
+{
+    rho_ = phases_[0]*phases_[0].thermo().rho();
+
+    for (label phasei=1; phasei<phases_.size(); phasei++)
+    {
+        rho_ += phases_[phasei]*phases_[phasei].thermo().rho();
+    }
+
+    forAll(phases_, phasei)
+    {
+        phases_[phasei].Alpha() =
+            phases_[phasei]*phases_[phasei].thermo().rho()/rho_;
+    }
+
+    calcAlphas();
+}
+
+
+void Foam::compressibleMultiphaseVoFMixture::correctRho
+(
+    const volScalarField& dp
+)
+{
+    forAll(phases_, phasei)
+    {
+        phases_[phasei].thermo().rho() += phases_[phasei].thermo().psi()*dp;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H
new file mode 100644
index 0000000000..5657d023fb
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H
@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::compressibleMultiphaseVoFMixture
+
+Description
+    Compressible multiphase mixture for interface-capturing simulations.
+
+    Derived from viscosity so that it can be used in conjunction with the
+    compressible turbulence models.
+
+    Surface tension and contact-angle is handled for the interface between each
+    phase-pair.
+
+SourceFiles
+    compressibleMultiphaseVoFMixture.C
+
+See also
+    Foam::multiphaseVoFMixture
+    Foam::compressibleVoFphase
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleMultiphaseVoFMixture_H
+#define compressibleMultiphaseVoFMixture_H
+
+#include "compressibleMultiphaseVoFMixtureThermo.H"
+#include "multiphaseVoFMixture.H"
+#include "viscosity.H"
+#include "UPtrListDictionary.H"
+#include "compressibleVoFphase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class compressibleMultiphaseVoFMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class compressibleMultiphaseVoFMixture
+:
+    public compressibleMultiphaseVoFMixtureThermo,
+    public multiphaseVoFMixture,
+    public viscosity
+{
+    // Private Data
+
+        //- Dictionary of phases
+        UPtrListDictionary<compressibleVoFphase> phases_;
+
+        volScalarField rho_;
+
+
+public:
+
+    TypeName("compressibleMultiphaseVoFMixture");
+
+
+    // Constructors
+
+        //- Construct from fvMesh
+        compressibleMultiphaseVoFMixture(const fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~compressibleMultiphaseVoFMixture()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        UPtrListDictionary<compressibleVoFphase>& phases()
+        {
+            return phases_;
+        }
+
+        //- Return the mixture density
+        virtual const volScalarField& rho() const
+        {
+            return rho_;
+        }
+
+        //- Return the kinematic laminar viscosity
+        virtual tmp<volScalarField> nu() const;
+
+        //- Return the laminar viscosity for patch
+        virtual tmp<scalarField> nu(const label patchi) const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        //- Correct the thermodynamics of each phase
+        virtual void correctThermo();
+
+        //- Return the effective temperature transport coefficient
+        //  derived from the phase internal energy equations i.e. <kappa/Cv>
+        virtual tmp<volScalarField> alphaEff
+        (
+            const volScalarField& nut
+        ) const;
+
+        //- Return the phase-averaged reciprocal Cv
+        tmp<volScalarField> rCv() const;
+
+        //- Update properties
+        virtual void correct();
+
+        //- Update densities for given pressure change
+        void correctRho(const volScalarField& dp);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/phase/phase.C b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixtureThermo/compressibleMultiphaseVoFMixtureThermo.C
similarity index 76%
rename from applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/phase/phase.C
rename to applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixtureThermo/compressibleMultiphaseVoFMixtureThermo.C
index 9a8bc613c7..26893351a9 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/phase/phase.C
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixtureThermo/compressibleMultiphaseVoFMixtureThermo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,17 +23,21 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "phase.H"
+#include "compressibleMultiphaseVoFMixtureThermo.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::phase::phase(const word& phaseName, const fvMesh& mesh)
+Foam::compressibleMultiphaseVoFMixtureThermo::
+compressibleMultiphaseVoFMixtureThermo
+(
+    const fvMesh& mesh
+)
 :
-    volScalarField
+    p_
     (
         IOobject
         (
-            IOobject::groupName("alpha", phaseName),
+            "p",
             mesh.time().name(),
             mesh,
             IOobject::MUST_READ,
@@ -41,7 +45,19 @@ Foam::phase::phase(const word& phaseName, const fvMesh& mesh)
         ),
         mesh
     ),
-    name_(phaseName)
+
+    T_
+    (
+        IOobject
+        (
+            "T",
+            mesh.time().name(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    )
 {}
 
 
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/phase/phase.H b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixtureThermo/compressibleMultiphaseVoFMixtureThermo.H
similarity index 70%
rename from applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/phase/phase.H
rename to applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixtureThermo/compressibleMultiphaseVoFMixtureThermo.H
index 5051e32ef1..4b5c1c6343 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/phase/phase.H
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixtureThermo/compressibleMultiphaseVoFMixtureThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,20 +22,17 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::phase
+    Foam::compressibleMultiphaseVoFMixtureThermo
 
 Description
-    Single incompressible phase derived from the phase-fraction.
-    Used as part of the multiPhaseMixture for interface-capturing multi-phase
-    simulations.
 
 SourceFiles
-    phase.C
+    compressibleMultiphaseVoFMixtureThermo.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef phase_H
-#define phase_H
+#ifndef compressibleMultiphaseVoFMixtureThermo_H
+#define compressibleMultiphaseVoFMixtureThermo_H
 
 #include "volFields.H"
 
@@ -45,40 +42,43 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                            Class phase Declaration
+                   Class compressibleMultiphaseVoFMixtureThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-class phase
-:
-    public volScalarField
+class compressibleMultiphaseVoFMixtureThermo
 {
+
+private:
+
     // Private Data
 
-        word name_;
+        //- Pressure
+        volScalarField p_;
+
+        //- Mixture temperature
+        volScalarField T_;
 
 
 public:
 
     // Constructors
 
-        //- Construct from components
-        phase
-        (
-            const word& name,
-            const fvMesh& mesh
-        );
+        //- Construct from fvMesh
+        compressibleMultiphaseVoFMixtureThermo(const fvMesh& mesh);
 
 
     // Member Functions
 
-        const word& name() const
+        //- Return pressure [Pa]
+        volScalarField& p()
         {
-            return name_;
+            return p_;
         }
 
-        const word& keyword() const
+        //- Return mixture temperature [K]
+        volScalarField& T()
         {
-            return name();
+            return T_;
         }
 };
 
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressiblePhase/compressiblePhase.C b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.C
similarity index 79%
rename from applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressiblePhase/compressiblePhase.C
rename to applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.C
index b7b7182280..0d3cf2721b 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressiblePhase/compressiblePhase.C
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,27 +23,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "compressiblePhase.H"
+#include "compressibleVoFphase.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::compressiblePhase::compressiblePhase
+Foam::compressibleVoFphase::compressibleVoFphase
 (
     const word& name,
+    const fvMesh& mesh,
     const volScalarField& T
 )
 :
-    phase(name, T.mesh()),
+    VoFphase(name, mesh),
     thermo_(nullptr),
     Alpha_
     (
         IOobject
         (
             IOobject::groupName("Alpha", name),
-            T.mesh().time().name(),
-            T.mesh()
+            mesh.time().name(),
+            mesh
         ),
-        T.mesh(),
+        mesh,
         dimensionedScalar(dimless, 0)
     ),
     dgdt_
@@ -51,12 +52,12 @@ Foam::compressiblePhase::compressiblePhase
         IOobject
         (
             IOobject::groupName("dgdt", name),
-            T.mesh().time().name(),
-            T.mesh(),
+            mesh.time().name(),
+            mesh,
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        T.mesh(),
+        mesh,
         dimensionedScalar(dimless/dimTime, 0)
     )
 {
@@ -65,21 +66,21 @@ Foam::compressiblePhase::compressiblePhase
         Tp.write();
     }
 
-    thermo_ = rhoThermo::New(T.mesh(), name);
+    thermo_ = rhoThermo::New(mesh, name);
     thermo_->validate(name, "e");
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::autoPtr<Foam::compressiblePhase> Foam::compressiblePhase::clone() const
+Foam::autoPtr<Foam::VoFphase> Foam::compressibleVoFphase::clone() const
 {
     NotImplemented;
-    return autoPtr<compressiblePhase>(nullptr);
+    return autoPtr<VoFphase>(nullptr);
 }
 
 
-void Foam::compressiblePhase::correct
+void Foam::compressibleVoFphase::correct
 (
     const volScalarField& p,
     const volScalarField& T
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressiblePhase/compressiblePhase.H b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.H
similarity index 71%
rename from applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressiblePhase/compressiblePhase.H
rename to applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.H
index 0499ef3e7c..bbc8fe982f 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressiblePhase/compressiblePhase.H
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleVoFphase/compressibleVoFphase.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,23 +22,27 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::compressiblePhase
+    Foam::compressibleVoFphase
 
 Description
-    Single incompressible phase derived from the phase-fraction.
-    Used as part of the multiPhaseMixture for interface-capturing multi-phase
-    simulations.
+    Single compressible phase derived from the VoFphase.
+
+    Used in compressibleMultiphaseVoFMixture for multiphase
+    interface-capturing simulations.
 
 SourceFiles
-    compressiblePhase.C
+    compressibleVoFphase.C
+
+See also
+    Foam::VoFphase
+    Foam::compressibleMultiphaseVoFMixture
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef compressiblePhase_H
-#define compressiblePhase_H
+#ifndef compressibleVoFphase_H
+#define compressibleVoFphase_H
 
-#include "phase.H"
-#include "dictionaryEntry.H"
+#include "VoFphase.H"
 #include "rhoThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -47,17 +51,22 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class compressiblePhase Declaration
+                            Class compressibleVoFphase Declaration
 \*---------------------------------------------------------------------------*/
 
-class compressiblePhase
+class compressibleVoFphase
 :
-    public phase
+    public VoFphase
 {
     // Private Data
 
+        //- Phase thermo
         autoPtr<rhoThermo> thermo_;
+
+        //- Phase mass-fraction
         volScalarField Alpha_;
+
+        //- Phase compressibility contribution
         volScalarField::Internal dgdt_;
 
 
@@ -66,36 +75,39 @@ public:
     // Constructors
 
         //- Construct from components
-        compressiblePhase
+        compressibleVoFphase
         (
-            const word& phaseName,
+            const word& name,
+            const fvMesh& mesh,
             const volScalarField& T
         );
 
         //- Return clone
-        autoPtr<compressiblePhase> clone() const;
+        virtual autoPtr<VoFphase> clone() const;
 
-        //- Return a pointer to a new compressiblePhase created on freestore
-        //  from Istream
+        //- Return a pointer to a new compressibleVoFphase
+        //  created on freestore from Istream
         class iNew
+        :
+            public VoFphase::iNew
         {
             const volScalarField& T_;
 
         public:
 
-            iNew
-            (
-                const volScalarField& T
-            )
+            iNew(const fvMesh& mesh, const volScalarField& T)
             :
+                VoFphase::iNew(mesh),
                 T_(T)
             {}
 
-            autoPtr<compressiblePhase> operator()(Istream& is) const
+            virtual autoPtr<VoFphase> operator()(Istream& is) const
             {
-                return autoPtr<compressiblePhase>
+                const word name(is);
+
+                return autoPtr<VoFphase>
                 (
-                    new compressiblePhase(is, T_)
+                    new compressibleVoFphase(name, mesh_, T_)
                 );
             }
         };
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/momentumPredictor.C b/applications/solvers/modules/compressibleMultiphaseVoF/momentumPredictor.C
new file mode 100644
index 0000000000..b4fbf4132a
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/momentumPredictor.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleMultiphaseVoF.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::solvers::compressibleMultiphaseVoF::divDevTau
+(
+    volVectorField& U
+)
+{
+    return momentumTransport.divDevTau(U) - fvm::Sp(contErr(), U);
+}
+
+
+void Foam::solvers::compressibleMultiphaseVoF::momentumPredictor()
+{
+    multiphaseVoFSolver::momentumPredictor();
+
+    if (pimple.momentumPredictor())
+    {
+        K = 0.5*magSqr(U);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/pressureCorrector.C b/applications/solvers/modules/compressibleMultiphaseVoF/pressureCorrector.C
new file mode 100644
index 0000000000..c0a7389d51
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/pressureCorrector.C
@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleMultiphaseVoF.H"
+#include "constrainHbyA.H"
+#include "constrainPressure.H"
+#include "adjustPhi.H"
+#include "findRefCell.H"
+#include "fvcMeshPhi.H"
+#include "fvcFlux.H"
+#include "fvcDdt.H"
+#include "fvcSnGrad.H"
+#include "fvcSup.H"
+#include "fvcReconstruct.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleMultiphaseVoF::pressureCorrector()
+{
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    if (rAU.valid())
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    const surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+
+    while (pimple.correct())
+    {
+        const volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
+            fvc::flux(HbyA)
+          + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
+        );
+
+        MRF.makeRelative(phiHbyA);
+
+        const surfaceScalarField phig
+        (
+            (
+                mixture.surfaceTensionForce(U)
+              - buoyancy.ghf*fvc::snGrad(rho)
+            )*rAUf*mesh.magSf()
+        );
+
+        phiHbyA += phig;
+
+        // Update the pressure BCs to ensure flux consistency
+        constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
+
+        PtrList<fvScalarMatrix> p_rghEqnComps(phases.size());
+
+        forAll(phases, phasei)
+        {
+            const compressibleVoFphase& phase = phases[phasei];
+            const rhoThermo& thermo = phase.thermo();
+            const volScalarField& rho = thermo.rho()();
+
+            p_rghEqnComps.set
+            (
+                phasei,
+                (
+                    fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh))
+                  + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho)
+                  - (fvModels().source(phase, rho)&rho)
+                ).ptr()
+            );
+        }
+
+        // Cache p_rgh prior to solve for density update
+        const volScalarField p_rgh_0(p_rgh);
+
+        while (pimple.correctNonOrthogonal())
+        {
+            fvScalarMatrix p_rghEqnIncomp
+            (
+                fvc::div(phiHbyA)
+              - fvm::laplacian(rAUf, p_rgh)
+            );
+
+            tmp<fvScalarMatrix> p_rghEqnComp;
+
+            forAll(phases, phasei)
+            {
+                const compressibleVoFphase& phase = phases[phasei];
+
+                tmp<fvScalarMatrix> p_rghEqnCompi
+                (
+                    (max(phase, scalar(0))/phase.thermo().rho())
+                   *p_rghEqnComps[phasei]
+                );
+
+                if (phasei == 0)
+                {
+                    p_rghEqnComp = p_rghEqnCompi;
+                }
+                else
+                {
+                    p_rghEqnComp.ref() += p_rghEqnCompi;
+                }
+            }
+
+            solve(p_rghEqnComp + p_rghEqnIncomp);
+
+            if (pimple.finalNonOrthogonalIter())
+            {
+                forAll(phases, phasei)
+                {
+                    compressibleVoFphase& phase = phases[phasei];
+
+                    phase.dgdt() =
+                        pos0(phase)
+                       *(p_rghEqnComps[phasei] & p_rgh)/phase.thermo().rho();
+                }
+
+                phi = phiHbyA + p_rghEqnIncomp.flux();
+
+                p = max(p_rgh + mixture.rho()*buoyancy.gh, pMin);
+                p_rgh = p - rho*buoyancy.gh;
+                p_rgh.correctBoundaryConditions();
+
+                U = HbyA
+                  + rAU()*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+                U.correctBoundaryConditions();
+                fvConstraints().constrain(U);
+            }
+        }
+
+        // Correct Uf if the mesh is moving
+        fvc::correctUf(Uf, U, fvc::absolute(phi, U), MRF);
+
+        // Update densities from change in p_rgh
+        mixture.correctRho(p_rgh - p_rgh_0);
+        mixture.correct();
+
+        // Correct p_rgh for consistency with p and the updated densities
+        p_rgh = p - rho*buoyancy.gh;
+        p_rgh.correctBoundaryConditions();
+    }
+
+    K = 0.5*magSqr(U);
+
+    tUEqn.clear();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/thermophysicalPredictor.C b/applications/solvers/modules/compressibleMultiphaseVoF/thermophysicalPredictor.C
new file mode 100644
index 0000000000..e9ea0fe9e7
--- /dev/null
+++ b/applications/solvers/modules/compressibleMultiphaseVoF/thermophysicalPredictor.C
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleMultiphaseVoF.H"
+#include "fvcMeshPhi.H"
+#include "fvcDdt.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::solvers::compressibleMultiphaseVoF::thermophysicalPredictor()
+{
+    volScalarField& T = mixture.T();
+
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(rho, T) + fvm::div(rhoPhi, T) - fvm::Sp(contErr(), T)
+      - fvm::laplacian(mixture.alphaEff(momentumTransport.nut()), T)
+      + (
+            fvc::div(fvc::absolute(phi, U), p)()() // - contErr()/rho*p
+          + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()()
+          - (U()&(fvModels().source(rho, U)&U)()) - contErr()*K
+        )*mixture.rCv()()
+     ==
+        fvModels().source(rho, T)
+    );
+
+    TEqn.relax();
+
+    fvConstraints().constrain(TEqn);
+
+    TEqn.solve();
+
+    fvConstraints().constrain(T);
+
+    mixture.correctThermo();
+    mixture.correct();
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleVoF/Make/files b/applications/solvers/modules/compressibleVoF/Make/files
index 9ed053db6b..90cf57439d 100644
--- a/applications/solvers/modules/compressibleVoF/Make/files
+++ b/applications/solvers/modules/compressibleVoF/Make/files
@@ -3,6 +3,7 @@ compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C
 compressibleInterPhaseThermophysicalTransportModel/compressibleInterPhaseThermophysicalTransportModel.C
 
 alphaSuSp.C
+momentumPredictor.C
 thermophysicalPredictor.C
 pressureCorrector.C
 compressibleVoF.C
diff --git a/applications/solvers/modules/compressibleVoF/compressibleVoF.C b/applications/solvers/modules/compressibleVoF/compressibleVoF.C
index 391202602c..2c8d2e47ca 100644
--- a/applications/solvers/modules/compressibleVoF/compressibleVoF.C
+++ b/applications/solvers/modules/compressibleVoF/compressibleVoF.C
@@ -179,17 +179,6 @@ void Foam::solvers::compressibleVoF::prePredictor()
 }
 
 
-void Foam::solvers::compressibleVoF::momentumPredictor()
-{
-    twoPhaseVoFSolver::momentumPredictor();
-
-    if (pimple.momentumPredictor())
-    {
-        K = 0.5*magSqr(U);
-    }
-}
-
-
 void Foam::solvers::compressibleVoF::postCorrector()
 {
     if (pimple.correctTransport())
diff --git a/applications/solvers/modules/compressibleVoF/compressibleVoF.H b/applications/solvers/modules/compressibleVoF/compressibleVoF.H
index 38acf9dbb0..a51e0a6a8b 100644
--- a/applications/solvers/modules/compressibleVoF/compressibleVoF.H
+++ b/applications/solvers/modules/compressibleVoF/compressibleVoF.H
@@ -62,9 +62,6 @@ See also
 #include "compressibleTwoPhaseVoFMixture.H"
 #include "compressibleInterPhaseTransportModel.H"
 #include "compressibleInterPhaseThermophysicalTransportModel.H"
-#include "buoyancy.H"
-#include "pressureReference.H"
-#include "fvmSup.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -167,12 +164,7 @@ protected:
         );
 
         //- Return the momentum equation stress term
-        virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U)
-        {
-            return
-                momentumTransport.divDevTau(U)
-              - fvm::Sp(contErr1() + contErr2(), U);
-        }
+        virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U);
 
 
 public:
diff --git a/applications/solvers/modules/compressibleVoF/momentumPredictor.C b/applications/solvers/modules/compressibleVoF/momentumPredictor.C
new file mode 100644
index 0000000000..e5a7468415
--- /dev/null
+++ b/applications/solvers/modules/compressibleVoF/momentumPredictor.C
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleVoF.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::fvVectorMatrix> Foam::solvers::compressibleVoF::divDevTau
+(
+    volVectorField& U
+)
+{
+    return
+        momentumTransport.divDevTau(U)
+      - fvm::Sp(contErr1() + contErr2(), U);
+}
+
+
+void Foam::solvers::compressibleVoF::momentumPredictor()
+{
+    twoPhaseVoFSolver::momentumPredictor();
+
+    if (pimple.momentumPredictor())
+    {
+        K = 0.5*magSqr(U);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/modules/compressibleVoF/pressureCorrector.C b/applications/solvers/modules/compressibleVoF/pressureCorrector.C
index cfab605fda..f7fcf82099 100644
--- a/applications/solvers/modules/compressibleVoF/pressureCorrector.C
+++ b/applications/solvers/modules/compressibleVoF/pressureCorrector.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,8 +31,9 @@ License
 #include "fvcFlux.H"
 #include "fvcDdt.H"
 #include "fvcSnGrad.H"
-#include "fvcReconstruct.H"
 #include "fvmDiv.H"
+#include "fvmSup.H"
+#include "fvcReconstruct.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
@@ -55,13 +56,13 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
         rAU = 1.0/UEqn.A();
     }
 
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+    const surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
 
     const surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
 
     while (pimple.correct())
     {
-        volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
+        const volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
         surfaceScalarField phiHbyA
         (
             "phiHbyA",
@@ -71,7 +72,7 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
 
         MRF.makeRelative(phiHbyA);
 
-        surfaceScalarField phig
+        const surfaceScalarField phig
         (
             (
                 interface.surfaceTensionForce()
@@ -178,7 +179,7 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
         }
 
         // Cache p_rgh prior to solve for density update
-        volScalarField p_rgh_0(p_rgh);
+        const volScalarField p_rgh_0(p_rgh);
 
         while (pimple.correctNonOrthogonal())
         {
diff --git a/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C b/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C
index 9dc98b197e..c3bf4d3b53 100644
--- a/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C
+++ b/applications/solvers/modules/compressibleVoF/thermophysicalPredictor.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ License
 #include "fvcMeshPhi.H"
 #include "fvcDdt.H"
 #include "fvmDiv.H"
+#include "fvmSup.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.C b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.C
index d76c24e968..58b8e7a717 100644
--- a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.C
+++ b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.C
@@ -34,6 +34,26 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::incompressibleMultiphaseVoFMixture::mu() const
+{
+    tmp<volScalarField> tmu
+    (
+        phases_[0]*phases_[0].rho()*phases_[0].nu()
+    );
+    volScalarField& mu = tmu.ref();
+
+    for (label phasei=1; phasei<phases_.size(); phasei++)
+    {
+        mu += phases_[phasei]*phases_[phasei].rho()*phases_[phasei].nu();
+    }
+
+    return tmu;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::incompressibleMultiphaseVoFMixture::incompressibleMultiphaseVoFMixture
@@ -78,68 +98,6 @@ Foam::incompressibleMultiphaseVoFMixture::incompressibleMultiphaseVoFMixture
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField>
-Foam::incompressibleMultiphaseVoFMixture::mu() const
-{
-    tmp<volScalarField> tmu
-    (
-        phases_[0]*phases_[0].rho()*phases_[0].nu()
-    );
-    volScalarField& mu = tmu.ref();
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        mu += phases_[phasei]*phases_[phasei].rho()*phases_[phasei].nu();
-    }
-
-    return tmu;
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::incompressibleMultiphaseVoFMixture::mu(const label patchi) const
-{
-    tmp<scalarField> tmu
-    (
-        phases_[0].boundaryField()[patchi]
-       *phases_[0].rho().value()
-       *phases_[0].nu(patchi)
-    );
-    scalarField& mu = tmu.ref();
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        mu +=
-            phases_[phasei].boundaryField()[patchi]
-           *phases_[phasei].rho().value()
-           *phases_[phasei].nu(patchi);
-    }
-
-    return tmu;
-}
-
-
-Foam::tmp<Foam::surfaceScalarField>
-Foam::incompressibleMultiphaseVoFMixture::muf() const
-{
-    tmp<surfaceScalarField> tmuf
-    (
-        fvc::interpolate(phases_[0])
-       *phases_[0].rho()*fvc::interpolate(phases_[0].nu())
-    );
-    surfaceScalarField& muf = tmuf.ref();
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        muf +=
-            fvc::interpolate(phases_[phasei])
-           *phases_[phasei].rho()*fvc::interpolate(phases_[phasei].nu());
-    }
-
-    return tmuf;
-}
-
-
 Foam::tmp<Foam::volScalarField>
 Foam::incompressibleMultiphaseVoFMixture::nu() const
 {
@@ -154,13 +112,6 @@ Foam::incompressibleMultiphaseVoFMixture::nu(const label patchi) const
 }
 
 
-Foam::tmp<Foam::surfaceScalarField>
-Foam::incompressibleMultiphaseVoFMixture::nuf() const
-{
-    return muf()/fvc::interpolate(rho());
-}
-
-
 void Foam::incompressibleMultiphaseVoFMixture::correct()
 {
     forAll(phases_, phasei)
@@ -176,7 +127,7 @@ void Foam::incompressibleMultiphaseVoFMixture::correct()
     }
 
     // Update the mixture kinematic viscosity
-    nu_ = mu()/rho();
+    nu_ = mu()/rho_;
 
     calcAlphas();
 }
diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.H b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.H
index d239e83eb3..e43df61dd6 100644
--- a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.H
+++ b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleMultiphaseVoFMixture.H
@@ -76,6 +76,12 @@ class incompressibleMultiphaseVoFMixture
         volScalarField nu_;
 
 
+    // Private Member Functions
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+
 public:
 
     TypeName("incompressibleMultiphaseVoFMixture");
@@ -106,24 +112,12 @@ public:
             return rho_;
         }
 
-        //- Return the dynamic laminar viscosity
-        tmp<volScalarField> mu() const;
-
-        //- Return the dynamic laminar viscosity for patch
-        tmp<scalarField> mu(const label patchi) const;
-
-        //- Return the face-interpolated dynamic laminar viscosity
-        tmp<surfaceScalarField> muf() const;
-
         //- Return the kinematic laminar viscosity
         virtual tmp<volScalarField> nu() const;
 
         //- Return the laminar viscosity for patch
         virtual tmp<scalarField> nu(const label patchi) const;
 
-        //- Return the face-interpolated dynamic laminar viscosity
-        tmp<surfaceScalarField> nuf() const;
-
         //- Correct the mixture properties
         virtual void correct();
 
diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.C b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.C
index 49a6886ddd..3bed49dfee 100644
--- a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.C
+++ b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.C
@@ -54,12 +54,4 @@ void Foam::incompressibleVoFphase::correct()
 }
 
 
-bool Foam::incompressibleVoFphase::read(const dictionary& dict)
-{
-    dict.lookup("rho") >> rho_;
-
-    return true;
-}
-
-
 // ************************************************************************* //
diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.H b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.H
index 9cad0da38b..31753ff583 100644
--- a/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.H
+++ b/applications/solvers/modules/incompressibleMultiphaseVoF/incompressibleMultiphaseVoFMixture/incompressibleVoFphase/incompressibleVoFphase.H
@@ -134,9 +134,6 @@ public:
 
         //- Correct the incompressibleVoFphase properties
         void correct();
-
-        //- Read base incompressibleVoFphaseProperties dictionary
-        bool read(const dictionary& incompressibleVoFphaseDict);
 };
 
 
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean
deleted file mode 100755
index a1028e7612..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean
+++ /dev/null
@@ -1,7 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-wclean libso compressibleMultiphaseMixture
-wclean
-
-#------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake
deleted file mode 100755
index 1cacde1c54..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake
+++ /dev/null
@@ -1,10 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # Run from this directory
-
-# Parse arguments for library compilation
-. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
-
-wmake $targetType compressibleMultiphaseMixture
-wmake $targetType
-
-#------------------------------------------------------------------------------
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files
deleted file mode 100644
index 8e6a4cb0de..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-compressibleMultiphaseInterFoam.C
-
-EXE = $(FOAM_APPBIN)/compressibleMultiphaseInterFoam
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
deleted file mode 100644
index 950197e1e9..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
+++ /dev/null
@@ -1,25 +0,0 @@
-EXE_INC = \
-    -IcompressibleMultiphaseMixture/lnInclude \
-    -I$(LIB_SRC)/physicalProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/VoF \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceCompression/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-    -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -lcompressibleMultiphaseMixture \
-    -lphysicalProperties \
-    -lfluidThermophysicalModels \
-    -lspecie \
-    -linterfaceCompression \
-    -linterfaceProperties \
-    -lmomentumTransportModels \
-    -lcompressibleMomentumTransportModels \
-    -lfiniteVolume \
-    -lfvModels \
-    -lfvConstraints \
-    -lmeshTools
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
deleted file mode 100644
index 30bddd01a7..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
+++ /dev/null
@@ -1,25 +0,0 @@
-{
-    fvScalarMatrix TEqn
-    (
-        fvm::ddt(rho, T) + fvm::div(mixture.rhoPhi(), T) - fvm::Sp(contErr, T)
-      - fvm::laplacian(mixture.alphaEff(turbulence->nut()), T)
-      + (
-            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
-          + (fvc::ddt(rho, K) + fvc::div(mixture.rhoPhi(), K))()()
-          - (U()&(fvModels.source(rho, U)&U)()) - contErr*K
-        )*mixture.rCv()()
-     ==
-        fvModels.source(rho, T)
-    );
-
-    TEqn.relax();
-
-    fvConstraints.constrain(TEqn);
-
-    TEqn.solve();
-
-    fvConstraints.constrain(T);
-
-    mixture.correctThermo();
-    mixture.correct();
-}
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
deleted file mode 100644
index 936820c8a4..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
+++ /dev/null
@@ -1,32 +0,0 @@
-fvVectorMatrix UEqn
-(
-    fvm::ddt(rho, U) + fvm::div(mixture.rhoPhi(), U) - fvm::Sp(contErr, U)
-  + turbulence->divDevTau(U)
- ==
-    fvModels.source(rho, U)
-);
-
-UEqn.relax();
-
-fvConstraints.constrain(UEqn);
-
-if (pimple.momentumPredictor())
-{
-    solve
-    (
-        UEqn
-     ==
-        fvc::reconstruct
-        (
-            (
-                mixture.surfaceTensionForce()
-              - ghf*fvc::snGrad(rho)
-              - fvc::snGrad(p_rgh)
-            ) * mesh.magSf()
-        )
-    );
-
-    fvConstraints.constrain(U);
-
-    K = 0.5*magSqr(U);
-}
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
deleted file mode 100644
index 0bbb978ea1..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
+++ /dev/null
@@ -1,122 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    compressibleMultiphaseInterFoam
-
-Description
-    Solver for n compressible, non-isothermal immiscible fluids using a VOF
-    (volume of fluid) phase-fraction based interface capturing approach.
-
-    The momentum and other fluid properties are of the "mixture" and a single
-    momentum equation is solved.
-
-    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "compressibleMultiphaseMixture.H"
-#include "compressibleMomentumTransportModels.H"
-#include "pimpleControl.H"
-#include "fvModels.H"
-#include "fvConstraints.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "postProcess.H"
-
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createTimeControls.H"
-    #include "createFields.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
-
-    volScalarField& p = mixture.p();
-    volScalarField& T = mixture.T();
-
-    turbulence->validate();
-
-    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (pimple.run(runTime))
-    {
-        #include "readTimeControls.H"
-        #include "CourantNo.H"
-        #include "alphaCourantNo.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.userTimeName() << nl << endl;
-
-        // --- Pressure-velocity PIMPLE corrector loop
-        while (pimple.loop())
-        {
-            fvModels.correct();
-
-            mixture.solve();
-
-            #include "contErr.H"
-
-            if (pimple.predictTransport())
-            {
-                turbulence->predict();
-            }
-
-            #include "UEqn.H"
-            #include "TEqn.H"
-
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                #include "pEqn.H"
-            }
-
-            if (pimple.correctTransport())
-            {
-                turbulence->correct();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/Make/files
deleted file mode 100644
index 2cd54fb34b..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/Make/files
+++ /dev/null
@@ -1,5 +0,0 @@
-phase/phase.C
-compressiblePhase/compressiblePhase.C
-compressibleMultiphaseMixture.C
-
-LIB = $(FOAM_LIBBIN)/libcompressibleMultiphaseMixture
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/Make/options
deleted file mode 100644
index 5dca2a4b02..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/Make/options
+++ /dev/null
@@ -1,13 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/physicalProperties/lnInclude \
-    -I$(LIB_SRC)/multiphaseModels/multiphaseProperties/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-LIB_LIBS = \
-    -lfluidThermophysicalModels \
-    -lmultiphaseProperties \
-    -lspecie \
-    -lfiniteVolume \
-    -lmeshTools
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressibleMultiphaseMixture.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressibleMultiphaseMixture.C
deleted file mode 100644
index 8f6b2ed90b..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressibleMultiphaseMixture.C
+++ /dev/null
@@ -1,648 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "compressibleMultiphaseMixture.H"
-#include "alphaContactAngleFvPatchScalarField.H"
-#include "correctContactAngle.H"
-#include "Time.H"
-#include "subCycle.H"
-#include "MULES.H"
-#include "fvcDiv.H"
-#include "fvcGrad.H"
-#include "fvcSnGrad.H"
-#include "fvcFlux.H"
-#include "fvcMeshPhi.H"
-#include "surfaceInterpolate.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(compressibleMultiphaseMixture, 0);
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::compressibleMultiphaseMixture::calcAlphas()
-{
-    scalar level = 0.0;
-    alphas_ == 0.0;
-
-    forAll(phases_, phasei)
-    {
-        alphas_ += level*phases_[phasei];
-        level += 1.0;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::compressibleMultiphaseMixture::compressibleMultiphaseMixture
-(
-    const volVectorField& U,
-    const surfaceScalarField& phi
-)
-:
-    IOdictionary
-    (
-        IOobject
-        (
-            "phaseProperties",
-            U.time().constant(),
-            U.db(),
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    ),
-
-    mesh_(U.mesh()),
-
-    p_
-    (
-        IOobject
-        (
-            "p",
-            U.mesh().time().name(),
-            U.mesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        U.mesh()
-    ),
-
-    T_
-    (
-        IOobject
-        (
-            "T",
-            U.mesh().time().name(),
-            U.mesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        U.mesh()
-    ),
-
-    phases_(lookup("phases"), compressiblePhase::iNew(T_)),
-
-    rho_
-    (
-        IOobject
-        (
-            "rho",
-            U.mesh().time().name(),
-            U.mesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        U.mesh(),
-        dimensionedScalar("rho", dimDensity, 0)
-    ),
-
-    U_(U),
-
-    phi_(phi),
-
-    rhoPhi_
-    (
-        IOobject
-        (
-            "rhoPhi",
-            mesh_.time().name(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar(dimMass/dimTime, 0)
-    ),
-
-    alphas_
-    (
-        IOobject
-        (
-            "alphas",
-            mesh_.time().name(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar(dimless, 0)
-    ),
-
-    sigmas_(lookup("sigmas")),
-    dimSigma_(1, 0, -2, 0, 0),
-    deltaN_
-    (
-        "deltaN",
-        1e-8/pow(average(mesh_.V()), 1.0/3.0)
-    )
-{
-    calcAlphas();
-    alphas_.write();
-    correct();
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-void Foam::compressibleMultiphaseMixture::correctThermo()
-{
-    forAll(phases_, phasei)
-    {
-        phases_[phasei].correct(p_, T_);
-    }
-}
-
-
-void Foam::compressibleMultiphaseMixture::correct()
-{
-    rho_ = phases_[0]*phases_[0].thermo().rho();
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        rho_ += phases_[phasei]*phases_[phasei].thermo().rho();
-    }
-
-    forAll(phases_, phasei)
-    {
-        phases_[phasei].Alpha() =
-            phases_[phasei]*phases_[phasei].thermo().rho()/rho_;
-    }
-}
-
-
-void Foam::compressibleMultiphaseMixture::correctRho(const volScalarField& dp)
-{
-    forAll(phases_, phasei)
-    {
-        phases_[phasei].thermo().rho() += phases_[phasei].thermo().psi()*dp;
-    }
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::compressibleMultiphaseMixture::nu() const
-{
-    volScalarField mu(phases_[0].Alpha()*phases_[0].thermo().mu());
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        mu += phases_[phasei].Alpha()*phases_[phasei].thermo().mu();
-    }
-
-    return mu/rho_;
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::compressibleMultiphaseMixture::nu
-(
-    const label patchi
-) const
-{
-    scalarField mu
-    (
-        phases_[0].Alpha().boundaryField()[patchi]
-       *phases_[0].thermo().mu(patchi)
-    );
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        mu +=
-            phases_[phasei].Alpha().boundaryField()[patchi]
-           *phases_[phasei].thermo().mu(patchi);
-    }
-
-    return mu/rho_.boundaryField()[patchi];
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::compressibleMultiphaseMixture::alphaEff
-(
-    const volScalarField& nut
-) const
-{
-    tmp<volScalarField> talphaEff
-    (
-        phases_[0]
-       *(
-           phases_[0].thermo().kappa()
-         + phases_[0].thermo().rho()*phases_[0].thermo().Cp()*nut
-        )/phases_[0].thermo().Cv()
-    );
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        talphaEff.ref() +=
-            phases_[phasei]
-           *(
-               phases_[phasei].thermo().kappa()
-             + phases_[phasei].thermo().rho()*phases_[phasei].thermo().Cp()*nut
-            )/phases_[phasei].thermo().Cv();
-    }
-
-    return talphaEff;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::compressibleMultiphaseMixture::rCv() const
-{
-    tmp<volScalarField> trCv(phases_[0]/phases_[0].thermo().Cv());
-
-    for (label phasei=1; phasei<phases_.size(); phasei++)
-    {
-        trCv.ref() += phases_[phasei]/phases_[phasei].thermo().Cv();
-    }
-
-    return trCv;
-}
-
-
-Foam::tmp<Foam::surfaceScalarField>
-Foam::compressibleMultiphaseMixture::surfaceTensionForce() const
-{
-    tmp<surfaceScalarField> tstf
-    (
-        surfaceScalarField::New
-        (
-            "surfaceTensionForce",
-            mesh_,
-            dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), 0)
-        )
-    );
-
-    surfaceScalarField& stf = tstf.ref();
-
-    forAll(phases_, phasei)
-    {
-        const compressiblePhase& alpha1 = phases_[phasei];
-
-        for (label phasej = phasei+1; phasej<phases_.size(); phasej++)
-        {
-            const compressiblePhase& alpha2 = phases_[phasej];
-
-            sigmaTable::const_iterator sigma =
-                sigmas_.find(interfacePair(alpha1, alpha2));
-
-            if (sigma == sigmas_.end())
-            {
-                FatalErrorInFunction
-                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
-                    << " in list of sigma values"
-                    << exit(FatalError);
-            }
-
-            stf += dimensionedScalar(dimSigma_, sigma())
-               *fvc::interpolate(K(alpha1, alpha2))*
-                (
-                    fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
-                  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
-                );
-        }
-    }
-
-    return tstf;
-}
-
-
-void Foam::compressibleMultiphaseMixture::solve()
-{
-    const Time& runTime = mesh_.time();
-
-    const dictionary& alphaControls = mesh_.solution().solverDict("alpha");
-    label nAlphaSubCycles(alphaControls.lookup<label>("nAlphaSubCycles"));
-    scalar cAlpha(alphaControls.lookup<scalar>("cAlpha"));
-
-    if (nAlphaSubCycles > 1)
-    {
-        surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
-        dimensionedScalar totalDeltaT = runTime.deltaT();
-
-        List<volScalarField*> alphaPtrs(phases_.size());
-        forAll(phases_, phasei)
-        {
-            alphaPtrs[phasei] = &phases_[phasei];
-        }
-
-        for
-        (
-            subCycle<volScalarField, subCycleFields> alphaSubCycle
-            (
-                alphaPtrs,
-                nAlphaSubCycles
-            );
-            !(++alphaSubCycle).end();
-        )
-        {
-            solveAlphas(cAlpha);
-            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
-        }
-
-        rhoPhi_ = rhoPhiSum;
-    }
-    else
-    {
-        solveAlphas(cAlpha);
-    }
-
-    correct();
-}
-
-
-Foam::tmp<Foam::surfaceVectorField> Foam::compressibleMultiphaseMixture::nHatfv
-(
-    const volScalarField& alpha1,
-    const volScalarField& alpha2
-) const
-{
-    /*
-    // Cell gradient of alpha
-    volVectorField gradAlpha =
-        alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
-
-    // Interpolated face-gradient of alpha
-    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
-    */
-
-    surfaceVectorField gradAlphaf
-    (
-        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
-      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
-    );
-
-    // Face unit interface normal
-    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
-}
-
-
-Foam::tmp<Foam::surfaceScalarField> Foam::compressibleMultiphaseMixture::nHatf
-(
-    const volScalarField& alpha1,
-    const volScalarField& alpha2
-) const
-{
-    // Face unit interface normal flux
-    return nHatfv(alpha1, alpha2) & mesh_.Sf();
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::compressibleMultiphaseMixture::K
-(
-    const phase& alpha1,
-    const phase& alpha2
-) const
-{
-    tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
-
-    correctContactAngle
-    (
-        alpha1,
-        alpha2,
-        U_.boundaryField(),
-        deltaN_,
-        tnHatfv.ref().boundaryFieldRef()
-    );
-
-    // Simple expression for curvature
-    return -fvc::div(tnHatfv & mesh_.Sf());
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::compressibleMultiphaseMixture::nearInterface() const
-{
-    tmp<volScalarField> tnearInt
-    (
-        volScalarField::New
-        (
-            "nearInterface",
-            mesh_,
-            dimensionedScalar(dimless, 0)
-        )
-    );
-
-    forAll(phases_, phasei)
-    {
-        tnearInt.ref() = max
-        (
-            tnearInt(),
-            pos0(phases_[phasei] - 0.01)*pos0(0.99 - phases_[phasei])
-        );
-    }
-
-    return tnearInt;
-}
-
-
-void Foam::compressibleMultiphaseMixture::solveAlphas
-(
-    const scalar cAlpha
-)
-{
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phirb,alpha)");
-
-    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
-    phic = min(cAlpha*phic, max(phic));
-
-    UPtrList<const volScalarField> alphas(phases_.size());
-    PtrList<surfaceScalarField> alphaPhis(phases_.size());
-
-    forAll(phases_, phasei)
-    {
-        const compressiblePhase& alpha = phases_[phasei];
-
-        alphas.set(phasei, &alpha);
-
-        alphaPhis.set
-        (
-            phasei,
-            new surfaceScalarField
-            (
-                phi_.name() + alpha.name(),
-                fvc::flux
-                (
-                    phi_,
-                    alpha,
-                    alphaScheme
-                )
-            )
-        );
-
-        surfaceScalarField& alphaPhi = alphaPhis[phasei];
-
-        forAll(phases_, phasej)
-        {
-            compressiblePhase& alpha2 = phases_[phasej];
-
-            if (&alpha2 == &alpha) continue;
-
-            surfaceScalarField phir(phic*nHatf(alpha, alpha2));
-
-            alphaPhi += fvc::flux
-            (
-                -fvc::flux(-phir, alpha2, alpharScheme),
-                alpha,
-                alpharScheme
-            );
-        }
-
-        // Limit alphaPhi for each phase
-        MULES::limit
-        (
-            1.0/mesh_.time().deltaT().value(),
-            geometricOneField(),
-            alpha,
-            phi_,
-            alphaPhi,
-            zeroField(),
-            zeroField(),
-            oneField(),
-            zeroField(),
-            false
-        );
-    }
-
-    MULES::limitSum(alphas, alphaPhis, phi_);
-
-    rhoPhi_ = dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), 0);
-
-    volScalarField sumAlpha
-    (
-        IOobject
-        (
-            "sumAlpha",
-            mesh_.time().name(),
-            mesh_
-        ),
-        mesh_,
-        dimensionedScalar(dimless, 0)
-    );
-
-    volScalarField divU(fvc::div(fvc::absolute(phi_, U_)));
-
-    forAll(phases_, phasei)
-    {
-        compressiblePhase& alpha = phases_[phasei];
-
-        surfaceScalarField& alphaPhi = alphaPhis[phasei];
-
-        volScalarField::Internal Sp
-        (
-            IOobject
-            (
-                "Sp",
-                mesh_.time().name(),
-                mesh_
-            ),
-            mesh_,
-            dimensionedScalar(alpha.dgdt().dimensions(), 0)
-        );
-
-        volScalarField::Internal Su
-        (
-            IOobject
-            (
-                "Su",
-                mesh_.time().name(),
-                mesh_
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            divU.v()*min(alpha.v(), scalar(1))
-        );
-
-        {
-            const scalarField& dgdt = alpha.dgdt();
-
-            forAll(dgdt, celli)
-            {
-                if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
-                {
-                    Sp[celli] += dgdt[celli]*alpha[celli];
-                    Su[celli] -= dgdt[celli]*alpha[celli];
-                }
-                else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
-                {
-                    Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
-                }
-            }
-        }
-
-        forAll(phases_, phasej)
-        {
-            const compressiblePhase& alpha2 = phases_[phasej];
-
-            if (&alpha2 == &alpha) continue;
-
-            const scalarField& dgdt2 = alpha2.dgdt();
-
-            forAll(dgdt2, celli)
-            {
-                if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
-                {
-                    Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
-                    Su[celli] += dgdt2[celli]*alpha[celli];
-                }
-                else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
-                {
-                    Sp[celli] += dgdt2[celli]*alpha2[celli];
-                }
-            }
-        }
-
-        MULES::explicitSolve
-        (
-            geometricOneField(),
-            alpha,
-            alphaPhi,
-            Sp,
-            Su
-        );
-
-        rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*alphaPhi;
-
-        Info<< alpha.name() << " volume fraction, min, max = "
-            << alpha.weightedAverage(mesh_.V()).value()
-            << ' ' << min(alpha).value()
-            << ' ' << max(alpha).value()
-            << endl;
-
-        sumAlpha += alpha;
-    }
-
-    Info<< "Phase-sum volume fraction, min, max = "
-        << sumAlpha.weightedAverage(mesh_.V()).value()
-        << ' ' << min(sumAlpha).value()
-        << ' ' << max(sumAlpha).value()
-        << endl;
-
-    calcAlphas();
-}
-
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressibleMultiphaseMixture.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressibleMultiphaseMixture.H
deleted file mode 100644
index d4ae286905..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseMixture/compressibleMultiphaseMixture.H
+++ /dev/null
@@ -1,295 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::compressibleMultiphaseMixture
-
-Description
-
-SourceFiles
-    compressibleMultiphaseMixture.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef compressibleMultiphaseMixture_H
-#define compressibleMultiphaseMixture_H
-
-#include "compressiblePhase.H"
-#include "PtrListDictionary.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                   Class compressibleMultiphaseMixture Declaration
-\*---------------------------------------------------------------------------*/
-
-class compressibleMultiphaseMixture
-:
-    public IOdictionary,
-    public viscosity
-{
-public:
-
-    class interfacePair
-    :
-        public Pair<word>
-    {
-    public:
-
-        class hash
-        :
-            public Hash<interfacePair>
-        {
-        public:
-
-            hash()
-            {}
-
-            label operator()(const interfacePair& key) const
-            {
-                return word::hash()(key.first()) + word::hash()(key.second());
-            }
-        };
-
-
-        // Constructors
-
-            interfacePair()
-            {}
-
-            interfacePair(const word& alpha1Name, const word& alpha2Name)
-            :
-                Pair<word>(alpha1Name, alpha2Name)
-            {}
-
-            interfacePair(const phase& alpha1, const phase& alpha2)
-            :
-                Pair<word>(alpha1.name(), alpha2.name())
-            {}
-
-
-        // Friend Operators
-
-            friend bool operator==
-            (
-                const interfacePair& a,
-                const interfacePair& b
-            )
-            {
-                return
-                (
-                    ((a.first() == b.first()) && (a.second() == b.second()))
-                 || ((a.first() == b.second()) && (a.second() == b.first()))
-                );
-            }
-
-            friend bool operator!=
-            (
-                const interfacePair& a,
-                const interfacePair& b
-            )
-            {
-                return (!(a == b));
-            }
-    };
-
-
-private:
-
-    // Private Data
-
-        const fvMesh& mesh_;
-
-        //- Pressure
-        volScalarField p_;
-
-        //- Mixture temperature
-        volScalarField T_;
-
-        //- Dictionary of phases
-        PtrListDictionary<compressiblePhase> phases_;
-
-        volScalarField rho_;
-
-        const volVectorField& U_;
-        const surfaceScalarField& phi_;
-
-        surfaceScalarField rhoPhi_;
-
-        volScalarField alphas_;
-
-        typedef HashTable<scalar, interfacePair, interfacePair::hash>
-            sigmaTable;
-
-        sigmaTable sigmas_;
-        dimensionSet dimSigma_;
-
-        //- Stabilisation for normalisation of the interface normal
-        const dimensionedScalar deltaN_;
-
-
-    // Private Member Functions
-
-        void calcAlphas();
-
-        void solveAlphas(const scalar cAlpha);
-
-        tmp<surfaceVectorField> nHatfv
-        (
-            const volScalarField& alpha1,
-            const volScalarField& alpha2
-        ) const;
-
-        tmp<surfaceScalarField> nHatf
-        (
-            const volScalarField& alpha1,
-            const volScalarField& alpha2
-        ) const;
-
-        tmp<volScalarField> K
-        (
-            const phase& alpha1,
-            const phase& alpha2
-        ) const;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("compressibleMultiphaseMixture");
-
-
-    // Constructors
-
-        //- Construct from components
-        compressibleMultiphaseMixture
-        (
-            const volVectorField& U,
-            const surfaceScalarField& phi
-        );
-
-
-    //- Destructor
-    virtual ~compressibleMultiphaseMixture()
-    {}
-
-
-    // Member Functions
-
-        //- Return the phases
-        const PtrListDictionary<compressiblePhase>& phases() const
-        {
-            return phases_;
-        }
-
-        //- Return non-const access to the phases
-        PtrListDictionary<compressiblePhase>& phases()
-        {
-            return phases_;
-        }
-
-        //- Return pressure [Pa]
-        volScalarField& p()
-        {
-            return p_;
-        }
-
-        //- Return mixture temperature [K]
-        volScalarField& T()
-        {
-            return T_;
-        }
-
-        //- Return mixture density [kg/m^3]
-        const volScalarField& rho() const
-        {
-            return rho_;
-        }
-
-        //- Return mixture velocity
-        const volVectorField& U() const
-        {
-            return U_;
-        }
-
-        //- Return the volumetric flux
-        const surfaceScalarField& phi() const
-        {
-            return phi_;
-        }
-
-        const surfaceScalarField& rhoPhi() const
-        {
-            return rhoPhi_;
-        }
-
-        //- Correct the thermodynamics of each phase
-        virtual void correctThermo();
-
-        //- Update properties
-        virtual void correct();
-
-        //- Update densities for given pressure change
-        void correctRho(const volScalarField& dp);
-
-        //- Kinematic viscosity of mixture [m^2/s]
-        virtual tmp<volScalarField> nu() const;
-
-        //- Kinematic viscosity of mixture for patch [m^2/s]
-        virtual tmp<scalarField> nu(const label patchi) const;
-
-        //- Return the effective temperature transport coefficient
-        //  derived from the phase internal energy equations i.e. <kappa/Cv>
-        virtual tmp<volScalarField> alphaEff
-        (
-            const volScalarField& nut
-        ) const;
-
-        //- Return the phase-averaged reciprocal Cv
-        tmp<volScalarField> rCv() const;
-
-        tmp<surfaceScalarField> surfaceTensionForce() const;
-
-        //- Indicator of the proximity of the interface
-        //  Field values are 1 near and 0 away for the interface.
-        tmp<volScalarField> nearInterface() const;
-
-        //- Solve for the mixture phase-fractions
-        void solve();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/contErr.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/contErr.H
deleted file mode 100644
index 04c1ac2115..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/contErr.H
+++ /dev/null
@@ -1,10 +0,0 @@
-volScalarField::Internal contErr
-(
-    (fvc::ddt(rho) + fvc::div(mixture.rhoPhi()))()
-);
-
-forAll(mixture.phases(), phasei)
-{
-    const volScalarField& rho = mixture.phases()[phasei].thermo().rho();
-    contErr -= (fvModels.source(mixture.phases()[phasei], rho)&rho);
-}
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
deleted file mode 100644
index 8aa8d690d0..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
+++ /dev/null
@@ -1,62 +0,0 @@
-Info<< "Reading field p_rgh\n" << endl;
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.name(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "createPhi.H"
-
-Info<< "Constructing compressibleMultiphaseMixture\n" << endl;
-compressibleMultiphaseMixture mixture(U, phi);
-
-const volScalarField& rho = mixture.rho();
-
-dimensionedScalar pMin("pMin", dimPressure, mixture);
-
-mesh.schemes().setFluxRequired(p_rgh.name());
-
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-// Construct compressible turbulence model
-autoPtr<compressible::momentumTransportModel> turbulence
-(
-    compressible::momentumTransportModel::New
-    (
-        rho,
-        U,
-        mixture.rhoPhi(),
-        mixture
-    )
-);
-
-Info<< "Creating field kinetic energy K\n" << endl;
-volScalarField K("K", 0.5*magSqr(U));
-
-#include "createFvModels.H"
-#include "createFvConstraints.H"
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
deleted file mode 100644
index 2b5d3d011a..0000000000
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ /dev/null
@@ -1,114 +0,0 @@
-{
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::flux(HbyA)
-      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
-    );
-
-    surfaceScalarField phig
-    (
-        (
-            mixture.surfaceTensionForce()
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
-    );
-
-    phiHbyA += phig;
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf);
-
-    PtrList<fvScalarMatrix> p_rghEqnComps(mixture.phases().size());
-
-    forAll(mixture.phases(), phasei)
-    {
-        const compressiblePhase& phase = mixture.phases()[phasei];
-        const rhoThermo& thermo = phase.thermo();
-        const volScalarField& rho = thermo.rho()();
-
-        p_rghEqnComps.set
-        (
-            phasei,
-            (
-                fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh))
-              + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho)
-              - (fvModels.source(phase, rho)&rho)
-            ).ptr()
-        );
-    }
-
-    // Cache p_rgh prior to solve for density update
-    volScalarField p_rgh_0(p_rgh);
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix p_rghEqnIncomp
-        (
-            fvc::div(phiHbyA)
-          - fvm::laplacian(rAUf, p_rgh)
-        );
-
-        tmp<fvScalarMatrix> p_rghEqnComp;
-
-        forAll(mixture.phases(), phasei)
-        {
-            const compressiblePhase& phase = mixture.phases()[phasei];
-
-            tmp<fvScalarMatrix> p_rghEqnCompi
-            (
-                (max(phase, scalar(0))/phase.thermo().rho())
-               *p_rghEqnComps[phasei]
-            );
-
-            if (phasei == 0)
-            {
-                p_rghEqnComp = p_rghEqnCompi;
-            }
-            else
-            {
-                p_rghEqnComp.ref() += p_rghEqnCompi;
-            }
-        }
-
-        solve(p_rghEqnComp + p_rghEqnIncomp);
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            p = max(p_rgh + mixture.rho()*gh, pMin);
-            p_rgh = p - rho*gh;
-
-            forAll(mixture.phases(), phasei)
-            {
-                compressiblePhase& phase = mixture.phases()[phasei];
-
-                phase.dgdt() =
-                    pos0(phase)
-                  *(p_rghEqnComps[phasei] & p_rgh)/phase.thermo().rho();
-            }
-
-            phi = phiHbyA + p_rghEqnIncomp.flux();
-
-            U = HbyA
-              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
-            U.correctBoundaryConditions();
-            fvConstraints.constrain(U);
-        }
-    }
-
-    // Update densities from change in p_rgh
-    mixture.correctRho(p_rgh - p_rgh_0);
-    mixture.correct();
-
-    // Correct p_rgh for consistency with p and the updated densities
-    p_rgh = p - rho*gh;
-    p_rgh.correctBoundaryConditions();
-
-    K = 0.5*magSqr(U);
-
-    Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
-}
diff --git a/bin/compressibleMultiphaseInterFoam b/bin/compressibleMultiphaseInterFoam
new file mode 100755
index 0000000000..bc7f30c28f
--- /dev/null
+++ b/bin/compressibleMultiphaseInterFoam
@@ -0,0 +1,46 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     | Website:  https://openfoam.org
+#   \\  /    A nd           | Copyright (C) 2023 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM.
+#
+#     OpenFOAM is free software: you can redistribute it and/or modify it
+#     under the terms of the GNU General Public License as published by
+#     the Free Software Foundation, either version 3 of the License, or
+#     (at your option) any later version.
+#
+#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#     for more details.
+#
+#     You should have received a copy of the GNU General Public License
+#     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Script
+#     compressibleMultiphaseInterFoam
+#
+# Description
+#     Script to inform the user that compressibleMultiphaseInterFoam has been
+#     superseded and replaced by the more general compressibleMultiphaseVoF
+#     solver module executed by the foamRun application.
+#
+#------------------------------------------------------------------------------
+
+cat <<EOF
+
+compressibleMultiphaseInterFoam has been superseded and replaced by the more general
+compressibleMultiphaseVoF solver module executed by the foamRun application:
+
+    foamRun -solver compressibleMultiphaseVoF
+
+EOF
+
+exec env foamRun -solver compressibleMultiphaseVoF "$@"
+
+#------------------------------------------------------------------------------
diff --git a/bin/multiphaseInterFoam b/bin/multiphaseInterFoam
index 54056a95e6..ebf71b2959 100755
--- a/bin/multiphaseInterFoam
+++ b/bin/multiphaseInterFoam
@@ -26,8 +26,8 @@
 #     multiphaseInterFoam
 #
 # Description
-#     Script to inform the user that multiphaseInterFoam has been superseded
-#     and replaced by the more general incompressibleMultiphaseVoF solver module
+#     Script to inform the user that multiphaseInterFoam has been superseded and
+#     replaced by the more general incompressibleMultiphaseVoF solver module
 #     executed by the foamRun application.
 #
 #------------------------------------------------------------------------------
diff --git a/etc/config.sh/bash_completion b/etc/config.sh/bash_completion
index 682a82f230..f8a00949b5 100644
--- a/etc/config.sh/bash_completion
+++ b/etc/config.sh/bash_completion
@@ -495,31 +495,6 @@ _combinePatchFaces_ ()
 }
 complete -o filenames -o nospace -F _combinePatchFaces_ combinePatchFaces
 
-_compressibleMultiphaseInterFoam_ ()
-{
-    local cur="${COMP_WORDS[COMP_CWORD]}"
-    local prev="${COMP_WORDS[COMP_CWORD-1]}"
-    local line=${COMP_LINE}
-    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
-
-    opts="-case -doc -fileHandler -help -hostRoots -libs -listFunctionObjects -listFvConstraints -listFvModels -listMomentumTransportModels -listScalarBCs -listSwitches -listVectorBCs -noFunctionObjects -parallel -postProcess -roots -srcDoc"
-    for o in $used ; do opts="${opts/$o/}" ; done
-    extra=""
-
-    [ "$COMP_CWORD" = 1 ] || \
-    case "$prev" in
-        -case)
-            opts="" ; extra="-d" ;;
-        -fileHandler)
-            opts="uncollated collated masterUncollated" ; extra="" ;;
-        -hostRoots|-libs|-roots)
-            opts="" ; extra="" ;;
-       *) ;;
-    esac
-    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
-}
-complete -o filenames -o nospace -F _compressibleMultiphaseInterFoam_ compressibleMultiphaseInterFoam
-
 _createBaffles_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"
@@ -4437,6 +4412,31 @@ _compressibleInterFoam_ ()
 }
 complete -o filenames -o nospace -F _compressibleInterFoam_ compressibleInterFoam
 
+_compressibleMultiphaseInterFoam_ ()
+{
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    local line=${COMP_LINE}
+    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
+
+    opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -parallel -roots -solver -srcDoc"
+    for o in $used ; do opts="${opts/$o/}" ; done
+    extra=""
+
+    [ "$COMP_CWORD" = 1 ] || \
+    case "$prev" in
+        -case)
+            opts="" ; extra="-d" ;;
+        -fileHandler)
+            opts="uncollated collated masterUncollated" ; extra="" ;;
+        -hostRoots|-libs|-roots|-solver)
+            opts="" ; extra="" ;;
+       *) ;;
+    esac
+    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
+}
+complete -o filenames -o nospace -F _compressibleMultiphaseInterFoam_ compressibleMultiphaseInterFoam
+
 _foamCleanCase_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"
@@ -4955,6 +4955,31 @@ _multiphaseEulerFoam_ ()
 }
 complete -o filenames -o nospace -F _multiphaseEulerFoam_ multiphaseEulerFoam
 
+_multiphaseInterFoam_ ()
+{
+    local cur="${COMP_WORDS[COMP_CWORD]}"
+    local prev="${COMP_WORDS[COMP_CWORD-1]}"
+    local line=${COMP_LINE}
+    local used=$(echo "$line" | grep -oE "\-[a-zA-Z]+ ")
+
+    opts="-case -doc -fileHandler -help -hostRoots -libs -noFunctionObjects -parallel -roots -solver -srcDoc"
+    for o in $used ; do opts="${opts/$o/}" ; done
+    extra=""
+
+    [ "$COMP_CWORD" = 1 ] || \
+    case "$prev" in
+        -case)
+            opts="" ; extra="-d" ;;
+        -fileHandler)
+            opts="uncollated collated masterUncollated" ; extra="" ;;
+        -hostRoots|-libs|-roots|-solver)
+            opts="" ; extra="" ;;
+       *) ;;
+    esac
+    COMPREPLY=( $(compgen -W "${opts}" $extra -- ${cur}) )
+}
+complete -o filenames -o nospace -F _multiphaseInterFoam_ multiphaseInterFoam
+
 _paraFoam_ ()
 {
     local cur="${COMP_WORDS[COMP_CWORD]}"
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/T b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/T
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/T
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/T
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/U b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/U
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/U
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/U
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.air.orig b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.air.orig
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.air.orig
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.air.orig
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.mercury.orig b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.mercury.orig
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.mercury.orig
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.mercury.orig
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.oil.orig b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.oil.orig
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.oil.orig
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.oil.orig
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.water.orig b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.water.orig
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/alpha.water.orig
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/alpha.water.orig
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/p
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/p
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p_rgh b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/p_rgh
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p_rgh
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/0/p_rgh
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/Allclean
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allclean
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/Allclean
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/Allrun
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/Allrun
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/Allrun
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/g
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/g
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/g
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/momentumTransport b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/momentumTransport
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/momentumTransport
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/momentumTransport
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/phaseProperties b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/phaseProperties
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/phaseProperties
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/phaseProperties
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.air b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.air
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.air
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.air
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.mercury b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.mercury
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.mercury
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.mercury
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.oil b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.oil
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.oil
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.oil
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.water b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.water
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/constant/physicalProperties.water
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/constant/physicalProperties.water
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/blockMeshDict b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/blockMeshDict
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/blockMeshDict
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/blockMeshDict
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/controlDict
similarity index 94%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/controlDict
index f876ed50ac..7dda736747 100644
--- a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/controlDict
+++ b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/controlDict
@@ -14,7 +14,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     compressibleMultiphaseInterFoam;
+application     foamRun;
+
+solver          compressibleMultiphaseVoF;
 
 startFrom       startTime;
 
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/decomposeParDict
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/decomposeParDict
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/decomposeParDict
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/fvSchemes
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/fvSchemes
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/fvSolution
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSolution
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/fvSolution
diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict b/tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/setFieldsDict
similarity index 100%
rename from tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict
rename to tutorials/modules/compressibleMultiphaseVoF/damBreak4phaseLaminar/system/setFieldsDict