diff --git a/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.C b/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.C
index d8a2e1986199505a05d0023fde6b43b2896ccbb9..9255de6ac5c12b7a890e975ccecd6eea36432fd1 100644
--- a/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.C
+++ b/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.C
@@ -35,7 +35,8 @@ Foam::polynomialThermo<EquationOfState, PolySize>::polynomialThermo(Istream& is)
     EquationOfState(is),
     Hf_(readScalar(is)),
     cpPolynomial_("cpPolynomial", is),
-    dhPolynomial_("dhPolynomial", cpPolynomial_.integrate())
+    dhPolynomial_("dhPolynomial", cpPolynomial_.integrate()),
+    sPolynomial_("sPolynomial", cpPolynomial_.integrateMinus1())
 {}
 
 
@@ -48,7 +49,11 @@ Foam::Ostream& Foam::operator<<
     const polynomialThermo<EquationOfState, PolySize>& pt
 )
 {
-    os  << static_cast<const EquationOfState&>(pt) << tab << pt.Hf_;
+    os  << static_cast<const EquationOfState&>(pt) << tab
+        << pt.Hf_ << tab
+        << pt.cpPolynomial_ << tab
+        << pt.dhPolynomial_ << tab
+        << pt.sPolynomial;
 
     os.check
     (
diff --git a/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.H b/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.H
index 903cf1d2dbe5be41af753f90cf3e906b0fe4fae2..52ffaabd8f48a8caf965e90635951811fa50eb59 100644
--- a/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.H
+++ b/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermo.H
@@ -107,6 +107,9 @@ class polynomialThermo
         //- Enthalpy - derived from cp
         typename Polynomial<PolySize>::intPolyType dhPolynomial_;
 
+        //- Entropy - derived from cp
+        Polynomial<PolySize> sPolynomial_;
+
 
     // Private member functions
 
@@ -116,7 +119,8 @@ class polynomialThermo
             const EquationOfState& pt,
             const scalar Hf,
             const Polynomial<PolySize>& cpPoly,
-            const typename Polynomial<PolySize>::intPolyType& hPoly
+            const typename Polynomial<PolySize>::intPolyType& hPoly,
+            const Polynomial<PolySize>& sPoly
         );
 
 
diff --git a/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermoI.H b/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermoI.H
index a34f9fda8015f3c90e1a04c5b3a44ccb69b95a55..847b4a799838a2429c39c995e8d98cc8f44644b2 100644
--- a/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermoI.H
+++ b/src/thermophysicalModels/specie/thermo/polynomial/polynomialThermoI.H
@@ -34,13 +34,15 @@ inline Foam::polynomialThermo<EquationOfState, PolySize>::polynomialThermo
     const EquationOfState& pt,
     const scalar Hf,
     const Polynomial<PolySize>& cpPoly,
-    const typename Polynomial<PolySize>::intPolyType& dhPoly
+    const typename Polynomial<PolySize>::intPolyType& dhPoly,
+    const Polynomial<PolySize>& sPoly
 )
 :
     EquationOfState(pt),
     Hf_(Hf),
     cpPolynomial_(cpPoly),
-    dhPolynomial_(dhPoly)
+    dhPolynomial_(dhPoly),
+    sPolynomial_(sPoly)
 {}
 
 
@@ -56,7 +58,8 @@ inline Foam::polynomialThermo<EquationOfState, PolySize>::polynomialThermo
     EquationOfState(name, pt),
     Hf_(pt.Hf_),
     cpPolynomial_(pt.cpPolynomial_),
-    dhPolynomial_(pt.dhPolynomial_)
+    dhPolynomial_(pt.dhPolynomial_),
+    sPolynomial_(pt.sPolynomial_)
 {}
 
 
@@ -106,9 +109,13 @@ inline Foam::scalar Foam::polynomialThermo<EquationOfState, PolySize>::s
     const scalar T
 ) const
 {
-    notImplemented("scalar polynomialThermo<EquationOfState>::s");
+    WarningIn("scalar polynomialThermo<EquationOfState>::s(const scalar)")
+        << "Standard entropy of formation currently set to zero"
+        << endl;
 
-    return 0.0;
+    scalar Sf = 0.0;
+
+    return (sPolynomial_.evaluate(T) + Sf)*this->W();
 }
 
 
@@ -130,6 +137,7 @@ inline void Foam::polynomialThermo<EquationOfState, PolySize>::operator+=
     Hf_ = molr1*Hf_ + molr2*pt.Hf_;
     cpPolynomial_ = molr1*cpPolynomial_ + molr2*pt.cpPolynomial_;
     dhPolynomial_ = molr1*dhPolynomial_ + molr2*pt.dhPolynomial_;
+    sPolynomial_ = molr1*sPolynomial_ + molr2*pt.sPolynomial_;
 }
 
 
@@ -149,6 +157,7 @@ inline void Foam::polynomialThermo<EquationOfState, PolySize>::operator-=
     Hf_ = molr1*Hf_ - molr2*pt.Hf_;
     cpPolynomial_ = molr1*cpPolynomial_ - molr2*pt.cpPolynomial_;
     dhPolynomial_ = molr1*dhPolynomial_ - molr2*pt.dhPolynomial_;
+    sPolynomial_ = molr1*sPolynomial_ - molr2*pt.sPolynomial_;
 }
 
 
@@ -174,7 +183,8 @@ inline Foam::polynomialThermo<EquationOfState, PolySize> Foam::operator+
         eofs,
         molr1*pt1.Hf_ + molr2*pt2.Hf_,
         molr1*pt1.cpPolynomial_ + molr2*pt2.cpPolynomial_,
-        molr1*pt1.dhPolynomial_ + molr2*pt2.dhPolynomial_
+        molr1*pt1.dhPolynomial_ + molr2*pt2.dhPolynomial_,
+        molr1*pt1.sPolynomial_ + molr2*pt2.sPolynomial_
     );
 }
 
@@ -199,7 +209,8 @@ inline Foam::polynomialThermo<EquationOfState, PolySize> Foam::operator-
         eofs,
         molr1*pt1.Hf_ - molr2*pt2.Hf_,
         molr1*pt1.cpPolynomial_ - molr2*pt2.cpPolynomial_,
-        molr1*pt1.dhPolynomial_ - molr2*pt2.dhPolynomial_
+        molr1*pt1.dhPolynomial_ - molr2*pt2.dhPolynomial_,
+        molr1*pt1.sPolynomial_ - molr2*pt2.sPolynomial_
     );
 }
 
@@ -216,7 +227,8 @@ inline Foam::polynomialThermo<EquationOfState, PolySize> Foam::operator*
         s*static_cast<const EquationOfState&>(pt),
         pt.Hf_,
         pt.cpPolynomial_,
-        pt.dhPolynomial_
+        pt.dhPolynomial_,
+        pt.sPolynomial_
     );
 }