diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H
index 84c283e7a7..e7601ee5fb 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pEqnComps.H
@@ -2,6 +2,7 @@ PtrList<fvScalarMatrix> pEqnComps(phases.size());
 
 {
     PtrList<volScalarField> dmdts(fluid.dmdts());
+    PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
 
     forAll(phases, phasei)
     {
@@ -71,5 +72,9 @@ PtrList<fvScalarMatrix> pEqnComps(phases.size());
         {
             pEqnComp -= dmdts[phasei]/rho;
         }
+        if (d2mdtdps.set(phasei))
+        {
+            pEqnComp -= correction(fvm::Sp(d2mdtdps[phasei]/rho, p_rgh));
+        }
     }
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
index 272386b579..23a851c79b 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
@@ -68,7 +68,11 @@ ThermalPhaseChangePhaseSystem
 :
     BasePhaseSystem(mesh),
     volatile_(this->template lookupOrDefault<word>("volatile", "none")),
-    dmdt0s_(this->phases().size())
+    dmdt0s_(this->phases().size()),
+    pressureImplicit_
+    (
+        this->template lookupOrDefault<Switch>("pressureImplicit", true)
+    )
 {
     this->generatePairsAndSubModels
     (
@@ -137,6 +141,28 @@ ThermalPhaseChangePhaseSystem
             )
         );
 
+        d2mdtdpfs_.insert
+        (
+            pair,
+            new volScalarField
+            (
+                IOobject
+                (
+                    IOobject::groupName
+                    (
+                        "thermalPhaseChange:d2mdtdpf",
+                        pair.name()
+                    ),
+                    this->mesh().time().timeName(),
+                    this->mesh(),
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                this->mesh(),
+                dimensionedScalar((dimDensity/dimTime)/dimPressure, 0)
+            )
+        );
+
         Tfs_.insert
         (
             pair,
@@ -242,6 +268,25 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::dmdts() const
 }
 
 
+template<class BasePhaseSystem>
+Foam::PtrList<Foam::volScalarField>
+Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::d2mdtdps() const
+{
+    PtrList<volScalarField> d2mdtdps(BasePhaseSystem::d2mdtdps());
+
+    forAllConstIter(phaseSystem::dmdtfTable, d2mdtdpfs_, d2mdtdpfIter)
+    {
+        const phasePair& pair = this->phasePairs_[d2mdtdpfIter.key()];
+        const volScalarField& d2mdtdpf = *d2mdtdpfIter();
+
+        addField(pair.phase1(), "d2mdtdp", d2mdtdpf, d2mdtdps);
+        addField(pair.phase2(), "d2mdtdp", - d2mdtdpf, d2mdtdps);
+    }
+
+    return d2mdtdps;
+}
+
+
 template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::
@@ -498,6 +543,30 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
                   + pos(dmdtfNew)*phase2.Y(volatile_);
             }
 
+            if (pressureImplicit_)
+            {
+                volScalarField& d2mdtdpf(*this->d2mdtdpfs_[pair]);
+
+                const dimensionedScalar dp(rootSmall*thermo1.p().average());
+
+                const volScalarField dTsatdp
+                (
+                    (
+                        saturationModelIter()->Tsat(thermo1.p() + dp/2)
+                      - saturationModelIter()->Tsat(thermo1.p() - dp/2)
+                    )/dp
+                );
+
+                d2mdtdpf = (H1 + H2)*dTsatdp/L;
+
+                if (volatile_ != "none")
+                {
+                    d2mdtdpf *=
+                        neg0(dmdtfNew)*phase1.Y(volatile_)
+                      + pos(dmdtfNew)*phase2.Y(volatile_);
+                }
+            }
+
             H1 = this->heatTransferModels_[pair].first()->K();
             H2 = this->heatTransferModels_[pair].second()->K();
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
index 94e805fd5b..a594b8c963 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -101,6 +101,9 @@ class ThermalPhaseChangePhaseSystem
         //- Mass transfer rates
         phaseSystem::dmdtfTable dmdtfs_;
 
+        //- Mass transfer linearisation coeffs
+        phaseSystem::dmdtfTable d2mdtdpfs_;
+
         //- Interface temperatures
         phaseSystem::dmdtfTable Tfs_;
 
@@ -114,6 +117,9 @@ class ThermalPhaseChangePhaseSystem
         //  function
         PtrList<volScalarField> dmdt0s_;
 
+        //- Switch to control whether or not mass transfer rates are linearised
+        //  in the pressure equation
+        Switch pressureImplicit_;
 
     // Private Member Functions
 
@@ -144,6 +150,9 @@ public:
         //- Return the mass transfer rates for each phase
         virtual PtrList<volScalarField> dmdts() const;
 
+        //- Return the mass transfer linearisation coeffs for each phase
+        virtual PtrList<volScalarField> d2mdtdps() const;
+
         //- Return the momentum transfer matrices for the cell-based algorithm
         virtual autoPtr<phaseSystem::momentumTransferTable> momentumTransfer();
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
index 7bc08e3f57..8d05a2f768 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@@ -671,6 +671,12 @@ Foam::PtrList<Foam::volScalarField> Foam::phaseSystem::dmdts() const
 }
 
 
+Foam::PtrList<Foam::volScalarField> Foam::phaseSystem::d2mdtdps() const
+{
+    return PtrList<volScalarField>(phaseModels_.size());
+}
+
+
 bool Foam::phaseSystem::incompressible() const
 {
     forAll(phaseModels_, phasei)
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index 80e662083b..8264d34b8c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@@ -544,6 +544,10 @@ public:
             //- Return the mass transfer rates for each phase
             virtual PtrList<volScalarField> dmdts() const;
 
+            //- Return the mass transfer pressure implicit coefficients
+            //  for each phase
+            virtual PtrList<volScalarField> d2mdtdps() const;
+
             //- Return incompressibility
             bool incompressible() const;