diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
index 10e8496206..f33843b48f 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
@@ -9,7 +9,7 @@
       - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
-    pZones.addEnthalpySource(thermo, hEqn);
+    pZones.addEnthalpySource(thermo, rho, hEqn);
 
     hEqn.relax();
 
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
index 699d029e08..e9083bc13b 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
@@ -73,6 +73,7 @@ Foam::thermalPorousZone::thermalPorousZone
 void Foam::thermalPorousZone::addEnthalpySource
 (
     const basicThermo& thermo,
+    const volScalarField& rho,
     fvScalarMatrix& hEqn
 ) const
 {
@@ -85,7 +86,6 @@ void Foam::thermalPorousZone::addEnthalpySource
     const scalarField& V = mesh().V();
     scalarField& hDiag = hEqn.diag();
     scalarField& hSource = hEqn.source();
-    const scalarField& rho = thermo.rho();
 
     scalarField hZone = thermo.h(scalarField(cells.size(), T_.value()), cells);
     scalar rate = 1e6;
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
index ab4c982818..afcc79b30c 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
@@ -138,6 +138,7 @@ public:
         void addEnthalpySource
         (
             const basicThermo& thermo,
+            const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;
 };
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C
index 640aebc618..ec68f13d5c 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.C
@@ -50,12 +50,13 @@ Foam::thermalPorousZones::thermalPorousZones
 void Foam::thermalPorousZones::addEnthalpySource
 (
     const basicThermo& thermo,
+    const volScalarField& rho,
     fvScalarMatrix& hEqn
 ) const
 {
     forAll(*this, i)
     {
-        operator[](i).addEnthalpySource(thermo, hEqn);
+        operator[](i).addEnthalpySource(thermo, rho, hEqn);
     }
 }
 
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H
index e53698b3bf..fbe913c74a 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZones.H
@@ -92,6 +92,7 @@ public:
         void addEnthalpySource
         (
             const basicThermo& thermo,
+            const volScalarField& rho,
             fvScalarMatrix& hEqn
         ) const;
 };