diff --git a/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/Make/options b/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/Make/options
index 682bbd15a2..686d3a5b14 100644
--- a/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/Make/options
+++ b/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/saturationModels/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhaseMixture/lnInclude \
@@ -11,6 +12,7 @@ EXE_INC = \
 LIB_LIBS = \
     -lphysicalProperties \
     -lfluidThermophysicalModels \
+    -lsaturationModels \
     -lcompressibleTwoPhaseMixture \
     -lcompressibleCavitationModels \
     -lfiniteVolume \
diff --git a/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C b/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
index a5bc25f4f6..2cc5c3984a 100644
--- a/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
+++ b/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
@@ -149,7 +149,8 @@ void Foam::fv::compressible::VoFCavitation::addSup
         const volScalarField::Internal vDot2P(pCoeff*mDot12P[1]);
 
         eqn +=
-            (vDot2P - vDot1P)*(cavitation_->pSat() - rho*gh)
+            vDot2P*cavitation_->pSat1() - vDot1P*cavitation_->pSat2()
+          - (vDot2P - vDot1P)*rho*gh
           - fvm::Sp(vDot2P - vDot1P, eqn.psi());
     }
 }
diff --git a/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C b/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C
index bb81de1125..37d2ea1bd2 100644
--- a/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C
+++ b/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C
@@ -56,7 +56,7 @@ Foam::compressible::cavitationModels::Kunz::Kunz
     Cc_("Cc", dimless, dict),
     Cv_("Cv", dimless, dict),
 
-    p0_("0", pSat().dimensions(), 0.0)
+    p0_("0", dimPressure, 0)
 {
     correct();
 }
@@ -81,11 +81,14 @@ Foam::compressible::cavitationModels::Kunz::mDotcvAlphal() const
         min(max(alphal(), scalar(0)), scalar(1))
     );
 
+    const volScalarField::Internal pSatv(this->pSatv());
+    const volScalarField::Internal pSatl(this->pSatl());
+
     return Pair<tmp<volScalarField::Internal>>
     (
         mcCoeff_*sqr(limitedAlphal)
-       *max(p - pSat(), p0_)/max(p - pSat(), 0.01*pSat()),
-       -mvCoeff_*min(p - pSat(), p0_)
+       *max(p - pSatv, p0_)/max(p - pSatv, 0.01*pSatv),
+       -mvCoeff_*min(p - pSatl, p0_)
     );
 }
 
@@ -107,11 +110,14 @@ Foam::compressible::cavitationModels::Kunz::mDotcvP() const
         min(max(alphal(), scalar(0)), scalar(1))
     );
 
+    const volScalarField::Internal pSatv(this->pSatv());
+    const volScalarField::Internal pSatl(this->pSatl());
+
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*sqr(limitedAlphal)*(1.0 - limitedAlphal)
-       *pos0(p - pSat())/max(p - pSat(), 0.01*pSat()),
-        (-mvCoeff_)*limitedAlphal*neg(p - pSat())
+        mcCoeff_*sqr(limitedAlphal)*(1 - limitedAlphal)
+       *pos0(p - pSatv)/max(p - pSatv, 0.01*pSatv),
+       -mvCoeff_*limitedAlphal*neg(p - pSatl)
     );
 }
 
diff --git a/src/twoPhaseModels/compressibleCavitation/Make/options b/src/twoPhaseModels/compressibleCavitation/Make/options
index 3a9b9f8008..21522675ff 100644
--- a/src/twoPhaseModels/compressibleCavitation/Make/options
+++ b/src/twoPhaseModels/compressibleCavitation/Make/options
@@ -1,19 +1,17 @@
 EXE_INC = \
     -I$(LIB_SRC)/physicalProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/saturationModels/lnInclude \
+    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
-    -lcompressibleTwoPhaseMixture \
-    -lfluidThermophysicalModels \
-    -lspecie \
     -lphysicalProperties \
-    -ltwoPhaseMixture \
-    -ltwoPhaseProperties \
-    -linterfaceProperties \
+    -lfluidThermophysicalModels \
+    -lsaturationModels \
+    -lcompressibleTwoPhaseMixture \
     -lfiniteVolume \
     -lmeshTools
diff --git a/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C b/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C
index ec21192796..2724b2dbee 100644
--- a/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C
+++ b/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C
@@ -56,7 +56,7 @@ Foam::compressible::cavitationModels::Merkle::Merkle
     Cc_("Cc", dimless, dict),
     Cv_("Cv", dimless, dict),
 
-    p0_("0", pSat().dimensions(), 0.0),
+    p0_("0", dimPressure, 0),
 
     mcCoeff_(Cc_/(0.5*sqr(UInf_)*tInf_))
 {
@@ -79,8 +79,8 @@ Foam::compressible::cavitationModels::Merkle::mDotcvAlphal() const
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*max(p - pSat(), p0_),
-       -mvCoeff_*min(p - pSat(), p0_)
+        mcCoeff_*max(p - pSatv(), p0_),
+       -mvCoeff_*min(p - pSatl(), p0_)
     );
 }
 
@@ -103,8 +103,8 @@ Foam::compressible::cavitationModels::Merkle::mDotcvP() const
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*(1.0 - limitedAlphal)*pos0(p - pSat()),
-        (-mvCoeff_)*limitedAlphal*neg(p - pSat())
+        mcCoeff_*(1 - limitedAlphal)*pos0(p - pSatv()),
+       -mvCoeff_*limitedAlphal*neg(p - pSatl())
     );
 }
 
diff --git a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C
index 8bc3965404..aa3f2df255 100644
--- a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C
+++ b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C
@@ -62,7 +62,7 @@ Foam::compressible::cavitationModels::SchnerrSauer::SchnerrSauer
     Cc_("Cc", dimless, dict),
     Cv_("Cv", dimless, dict),
 
-    p0_("0", pSat().dimensions(), 0.0)
+    p0_("0", dimPressure, 0.0)
 {
     correct();
 }
@@ -79,7 +79,7 @@ Foam::compressible::cavitationModels::SchnerrSauer::rRb
     return pow
     (
         ((4*constant::mathematical::pi*n_)/3)
-       *limitedAlphal/(1.0 + alphaNuc() - limitedAlphal),
+       *limitedAlphal/(1 + alphaNuc() - limitedAlphal),
         1.0/3.0
     );
 }
@@ -96,7 +96,8 @@ Foam::compressible::cavitationModels::SchnerrSauer::alphaNuc() const
 Foam::tmp<Foam::volScalarField::Internal>
 Foam::compressible::cavitationModels::SchnerrSauer::pCoeff
 (
-    const volScalarField::Internal& p
+    const volScalarField::Internal& p,
+    const volScalarField::Internal& pSat
 ) const
 {
     const volScalarField::Internal limitedAlphal
@@ -106,13 +107,12 @@ Foam::compressible::cavitationModels::SchnerrSauer::pCoeff
 
     const volScalarField::Internal rho
     (
-        limitedAlphal*rhol()
-      + (scalar(1) - limitedAlphal)*rhov()
+        limitedAlphal*rhol() + (1 - limitedAlphal)*rhov()
     );
 
     return
-        (3*rhol()*rhov())*sqrt(2/(3*rhol()))
-       *rRb(limitedAlphal)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
+        (3*rhol()*rhov())*sqrt((2.0/3.0)/rhol())
+       *rRb(limitedAlphal)/(rho*sqrt(mag(p - pSat) + 0.01*pSat));
 }
 
 
@@ -122,17 +122,21 @@ Foam::compressible::cavitationModels::SchnerrSauer::mDotcvAlphal() const
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal pCoeff(this->pCoeff(p));
-
     const volScalarField::Internal limitedAlphal
     (
         min(max(alphal(), scalar(0)), scalar(1))
     );
 
+    const volScalarField::Internal pSatv(this->pSatv());
+    const volScalarField::Internal pSatl(this->pSatl());
+
     return Pair<tmp<volScalarField::Internal>>
     (
-        Cc_*limitedAlphal*pCoeff*max(p - pSat(), p0_),
-       -Cv_*(1.0 + alphaNuc() - limitedAlphal)*pCoeff*min(p - pSat(), p0_)
+        Cc_*limitedAlphal*pCoeff(p, pSatv)*max(p - pSatv, p0_),
+       -Cv_
+       *(1 + alphaNuc() - limitedAlphal)
+       *pCoeff(p, pSatl)
+       *min(p - pSatl, p0_)
     );
 }
 
@@ -143,19 +147,22 @@ Foam::compressible::cavitationModels::SchnerrSauer::mDotcvP() const
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal pCoeff(this->pCoeff(p));
-
     const volScalarField::Internal limitedAlphal
     (
         min(max(alphal(), scalar(0)), scalar(1))
     );
 
-    const volScalarField::Internal apCoeff(limitedAlphal*pCoeff);
+    const volScalarField::Internal pSatv(this->pSatv());
+    const volScalarField::Internal pSatl(this->pSatl());
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        Cc_*(1.0 - limitedAlphal)*pos0(p - pSat())*apCoeff,
-        (-Cv_)*(1.0 + alphaNuc() - limitedAlphal)*neg(p - pSat())*apCoeff
+        Cc_*(1 - limitedAlphal)*pos0(p - pSatv)*limitedAlphal*pCoeff(p, pSatv),
+       -Cv_
+       *(1 + alphaNuc() - limitedAlphal)
+       *neg(p - pSatl)
+       *limitedAlphal
+       *pCoeff(p, pSatl)
     );
 }
 
diff --git a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H
index 96ff88a761..4822c16a26 100644
--- a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H
+++ b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H
@@ -90,7 +90,8 @@ class SchnerrSauer
         //- Part of the condensation and vapourisation rates
         tmp<volScalarField::Internal> pCoeff
         (
-            const volScalarField::Internal& p
+            const volScalarField::Internal& p,
+            const volScalarField::Internal& pSat
         ) const;
 
 
diff --git a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C
index e941a22df2..c8182811b2 100644
--- a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C
+++ b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C
@@ -48,13 +48,16 @@ Foam::compressible::cavitationModel::cavitationModel
 :
     phases_(phases),
     liquidIndex_(phases.index(dict.lookup<word>("liquid"))),
-    pSat_("pSat", dimPressure, dict.lookup("pSat"))
+    saturationModel_
+    (
+        saturationPressureModel::New(dict.subDict("saturationPressure"))
+    )
 {}
 
 
 bool Foam::compressible::cavitationModel::read(const dictionary& dict)
 {
-    dict.lookup("pSat") >> pSat_;
+    saturationModel_.reset(saturationPressureModel::New(dict).ptr());
 
     return true;
 }
diff --git a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H
index 19b3181b29..d5f809d896 100644
--- a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H
+++ b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H
@@ -38,6 +38,7 @@ SourceFiles
 #include "compressibleTwoPhases.H"
 #include "fvMatricesFwd.H"
 #include "Pair.H"
+#include "saturationPressureModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,8 +63,8 @@ protected:
         //- Index of the liquid
         const bool liquidIndex_;
 
-        //- Saturation vapour pressure
-        dimensionedScalar pSat_;
+        //- The saturation pressure model
+        autoPtr<saturationPressureModel> saturationModel_;
 
 
     // Protected Member Functions
@@ -92,6 +93,30 @@ protected:
             return phases_.rho(!liquidIndex_);
         }
 
+        //- Return the liquid thermo
+        inline const rhoThermo& thermol() const
+        {
+            return phases_.thermo(liquidIndex_);
+        }
+
+        //- Return the vapour thermo
+        inline const rhoThermo& thermov() const
+        {
+            return phases_.thermo(!liquidIndex_);
+        }
+
+        //- Return the saturation vapour pressure for the liquid
+        inline tmp<volScalarField::Internal> pSatl() const
+        {
+            return saturationModel_->pSat(thermol().T()());
+        }
+
+        //- Return the saturation vapour pressure for the vapour
+        inline tmp<volScalarField::Internal> pSatv() const
+        {
+            return saturationModel_->pSat(thermov().T()());
+        }
+
 
 public:
 
@@ -139,10 +164,16 @@ public:
 
     // Member Functions
 
-        //- Return the saturation vapour pressure
-        inline const dimensionedScalar& pSat() const
+        //- Return the saturation vapour pressure for phase 1
+        inline tmp<volScalarField::Internal> pSat1() const
         {
-            return pSat_;
+            return liquidIndex_ ? pSatv() : pSatl();
+        }
+
+        //- Return the saturation vapour pressure for phase 2
+        inline tmp<volScalarField::Internal> pSat2() const
+        {
+            return liquidIndex_ ? pSatl() : pSatv();
         }
 
         //- Return the mass condensation and vaporisation rates as a
diff --git a/src/twoPhaseModels/incompressibleCavitation/Make/options b/src/twoPhaseModels/incompressibleCavitation/Make/options
index f51d79235b..7001bab5bc 100644
--- a/src/twoPhaseModels/incompressibleCavitation/Make/options
+++ b/src/twoPhaseModels/incompressibleCavitation/Make/options
@@ -3,10 +3,11 @@ EXE_INC = \
     -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/twoPhaseModels/incompressibleTwoPhaseMixture/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
+    -lphysicalProperties \
     -lincompressibleTwoPhaseMixture \
-    -lmeshTools \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lmeshTools