diff --git a/include/cantera/base/Solution.h b/include/cantera/base/Solution.h
index 91b79bc279..e1876e1510 100644
--- a/include/cantera/base/Solution.h
+++ b/include/cantera/base/Solution.h
@@ -106,11 +106,14 @@ class Solution : public std::enable_shared_from_this<Solution>
 
     //! Get the number of adjacent phases
     size_t nAdjacent() const {
-         return m_adjacent.size();
+        return m_adjacent.size();
     }
 
     //! Get the name of an adjacent phase by index
     string adjacentName(size_t i) const {
+        if (i < 0 || i >= m_adjacent.size()) {
+            throw CanteraError("Solution::adjacentName", "Invalid index {}.", i);
+        }
         return m_adjacent.at(i)->name();
     }
 
diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h
index 1234479470..1e7539ccf5 100644
--- a/include/cantera/thermo/Phase.h
+++ b/include/cantera/thermo/Phase.h
@@ -875,7 +875,7 @@ class Phase
     bool m_caseSensitiveSpecies = false;
 
     //! Vector of size m_kk, used as a temporary holding area.
-    mutable vector<double> m_tmpV;
+    mutable vector<double> m_workS;
 
 private:
     //! Find lowercase species name in m_speciesIndices when case sensitive
diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h
index c89e7aa5b8..5845f7aea2 100644
--- a/include/cantera/thermo/ThermoPhase.h
+++ b/include/cantera/thermo/ThermoPhase.h
@@ -527,15 +527,13 @@ class ThermoPhase : public Phase
     /**
      * Returns the amount of enthalpy per mole,
      * @f[
-     * \bar{h}(T, P, X_k) = \sum_k X_k \bar{h}_k(T)
+     * \hat{h} = \sum_k X_k \hat{h}_k
      * @f]
-     * The formula is written in terms of the partial molar enthalpies.
-     * @f$ \bar{h}_k(T, p, m_k) @f$.
-     * @relates getPartialMolarEnthalpies()
+     * @see getPartialMolarEnthalpies()
      */
     virtual double enthalpy_mole() const {
-        getPartialMolarEnthalpies(m_tmpV.data());
-        return mean_X(m_tmpV);
+        getPartialMolarEnthalpies(m_workS.data());
+        return mean_X(m_workS);
     }
 
     //! Molar internal energy. Units: J/kmol.
@@ -547,40 +545,38 @@ class ThermoPhase : public Phase
     /**
      * Returns the amount of entropy per mole,
      * @f[
-     * \bar{s}(T, P, X_k) = \sum_k X_k \bar{s}_k(T)
+     * \hat{s} = \sum_k X_k \hat{s}_k
      * @f]
-     * The formula is written in terms of the partial molar entropies.
-     * @f$ \bar{s}_k(T, p, m_k) @f$.
-     * @relates getPartialMolarEnthalpies()
+     * @see getPartialMolarEnthalpies()
      */
     virtual double entropy_mole() const {
-        getPartialMolarEntropies(m_tmpV.data());
-        return mean_X(m_tmpV);
+        getPartialMolarEntropies(m_workS.data());
+        return mean_X(m_workS);
     }
 
     //! Molar Gibbs function. Units: J/kmol.
     /*!
      * Returns the Gibbs free energy per mole,
      * @f[
-     * \bar{g}(T, P, X_k) = \sum_k X_k \mu_k(T)
+     * \hat{g} = \sum_k X_k \mu_k
      * @f]
-     * @relates getChemPotentials()
+     * @see getChemPotentials()
      */
     virtual double gibbs_mole() const {
-        getChemPotentials(m_tmpV.data());
-        return mean_X(m_tmpV);
+        getChemPotentials(m_workS.data());
+        return mean_X(m_workS);
     }
 
     //! Molar heat capacity at constant pressure. Units: J/kmol/K.
     /*!
      * @f[
-     * \bar{c}_p(T, P, X_k) = \sum_k X_k \bar{c}_{p,k}(T)
+     * \hat{c}_p = \sum_k X_k \hat{c}_{p,k}
      * @f]
-     * @relates getPartialMolarCp()
+     * @see getPartialMolarCp()
      */
     virtual double cp_mole() const {
-        getPartialMolarCp(m_tmpV.data());
-        return mean_X(m_tmpV);
+        getPartialMolarCp(m_workS.data());
+        return mean_X(m_workS);
     }
 
     //! Molar heat capacity at constant volume. Units: J/kmol/K.
diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp
index d37c6efc68..b00198b3c1 100644
--- a/src/thermo/HMWSoln.cpp
+++ b/src/thermo/HMWSoln.cpp
@@ -53,8 +53,8 @@ HMWSoln::HMWSoln(const string& inputFile, const string& id_) :
 
 double HMWSoln::relative_enthalpy() const
 {
-    getPartialMolarEnthalpies(m_tmpV.data());
-    double hbar = mean_X(m_tmpV);
+    getPartialMolarEnthalpies(m_workS.data());
+    double hbar = mean_X(m_workS);
     getEnthalpy_RT(m_gamma_tmp.data());
     for (size_t k = 0; k < m_kk; k++) {
         m_gamma_tmp[k] *= RT();
@@ -66,20 +66,20 @@ double HMWSoln::relative_enthalpy() const
 double HMWSoln::relative_molal_enthalpy() const
 {
     double L = relative_enthalpy();
-    getMoleFractions(m_tmpV.data());
+    getMoleFractions(m_workS.data());
     double xanion = 0.0;
     size_t kcation = npos;
     double xcation = 0.0;
     size_t kanion = npos;
     for (size_t k = 0; k < m_kk; k++) {
         if (charge(k) > 0.0) {
-            if (m_tmpV[k] > xanion) {
-                xanion = m_tmpV[k];
+            if (m_workS[k] > xanion) {
+                xanion = m_workS[k];
                 kanion = k;
             }
         } else if (charge(k) < 0.0) {
-            if (m_tmpV[k] > xcation) {
-                xcation = m_tmpV[k];
+            if (m_workS[k] > xcation) {
+                xcation = m_workS[k];
                 kcation = k;
             }
         }
@@ -144,8 +144,8 @@ void HMWSoln::getActivityConcentrations(double* c) const
 
 double HMWSoln::standardConcentration(size_t k) const
 {
-    getStandardVolumes(m_tmpV.data());
-    double mvSolvent = m_tmpV[0];
+    getStandardVolumes(m_workS.data());
+    double mvSolvent = m_workS[0];
     if (k > 0) {
         return m_Mnaught / mvSolvent;
     }
@@ -1104,7 +1104,7 @@ double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
 
 void HMWSoln::initLengths()
 {
-    m_tmpV.resize(m_kk, 0.0);
+    m_workS.resize(m_kk, 0.0);
     m_molalitiesCropped.resize(m_kk, 0.0);
 
     size_t maxCounterIJlen = 1 + (m_kk-1) * (m_kk-2) / 2;
@@ -3936,7 +3936,7 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
 void HMWSoln::printCoeffs() const
 {
     calcMolalities();
-    vector<double>& moleF = m_tmpV;
+    vector<double>& moleF = m_workS;
 
     // Update the coefficients wrt Temperature. Calculate the derivatives as well
     s_updatePitzer_CoeffWRTemp(2);
diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp
index 59623f930a..d4a90092d8 100644
--- a/src/thermo/IdealMolalSoln.cpp
+++ b/src/thermo/IdealMolalSoln.cpp
@@ -47,8 +47,8 @@ IdealMolalSoln::IdealMolalSoln(const string& inputFile, const string& id_) :
 
 double IdealMolalSoln::intEnergy_mole() const
 {
-    getPartialMolarIntEnergies(m_tmpV.data());
-    return mean_X(m_tmpV);
+    getPartialMolarIntEnergies(m_workS.data());
+    return mean_X(m_workS);
 }
 
 // ------- Mechanical Equation of State Properties ------------------------
diff --git a/src/thermo/PengRobinson.cpp b/src/thermo/PengRobinson.cpp
index 4bb7dd2139..82d0db4be6 100644
--- a/src/thermo/PengRobinson.cpp
+++ b/src/thermo/PengRobinson.cpp
@@ -129,8 +129,8 @@ double PengRobinson::pressure() const
 
 double PengRobinson::standardConcentration(size_t k) const
 {
-    getStandardVolumes(m_tmpV.data());
-    return 1.0 / m_tmpV[k];
+    getStandardVolumes(m_workS.data());
+    return 1.0 / m_workS[k];
 }
 
 void PengRobinson::getActivityCoefficients(double* ac) const
@@ -256,9 +256,9 @@ void PengRobinson::getPartialMolarEntropies(double* sbar) const
     // Using the identity : (hk - T*sk) = gk
     double T = temperature();
     getPartialMolarEnthalpies(sbar);
-    getChemPotentials(m_tmpV.data());
+    getChemPotentials(m_workS.data());
     for (size_t k = 0; k < m_kk; k++) {
-        sbar[k] = (sbar[k] - m_tmpV[k])/T;
+        sbar[k] = (sbar[k] - m_workS[k])/T;
     }
 }
 
@@ -267,9 +267,9 @@ void PengRobinson::getPartialMolarIntEnergies(double* ubar) const
     // u_i = h_i - p*v_i
     double p = pressure();
     getPartialMolarEnthalpies(ubar);
-    getPartialMolarVolumes(m_tmpV.data());
+    getPartialMolarVolumes(m_workS.data());
     for (size_t k = 0; k < m_kk; k++) {
-        ubar[k] = ubar[k] - p*m_tmpV[k];
+        ubar[k] = ubar[k] - p*m_workS[k];
     }
 }
 
diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp
index 090a47919d..5837ab7791 100644
--- a/src/thermo/Phase.cpp
+++ b/src/thermo/Phase.cpp
@@ -806,7 +806,7 @@ bool Phase::addSpecies(shared_ptr<Species> spec)
         m_y.push_back(0.0);
         m_ym.push_back(0.0);
     }
-    m_tmpV.push_back(0.0);
+    m_workS.push_back(0.0);
     invalidateCache();
     return true;
 }
diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp
index c3d53b13d5..b3e05de7c2 100644
--- a/src/thermo/RedlichKwongMFTP.cpp
+++ b/src/thermo/RedlichKwongMFTP.cpp
@@ -140,8 +140,8 @@ double RedlichKwongMFTP::pressure() const
 
 double RedlichKwongMFTP::standardConcentration(size_t k) const
 {
-    getStandardVolumes(m_tmpV.data());
-    return 1.0 / m_tmpV[k];
+    getStandardVolumes(m_workS.data());
+    return 1.0 / m_workS[k];
 }
 
 void RedlichKwongMFTP::getActivityCoefficients(double* ac) const
@@ -237,10 +237,10 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(double* hbar) const
     double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
 
     for (size_t k = 0; k < m_kk; k++) {
-        m_tmpV[k] = 0.0;
+        m_workS[k] = 0.0;
         for (size_t i = 0; i < m_kk; i++) {
             size_t counter = k + m_kk*i;
-            m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
+            m_workS[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
         }
     }
 
@@ -248,7 +248,7 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(double* hbar) const
     double fac2 = mv + TKelvin * dpdT_ / dpdV_;
     for (size_t k = 0; k < m_kk; k++) {
         double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
-                       + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
+                       + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_workS[k]
                        +  b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
         hbar[k] = hbar[k] + hE_v;
         hbar[k] -= fac2 * dpdni_[k];
@@ -276,10 +276,10 @@ void RedlichKwongMFTP::getPartialMolarEntropies(double* sbar) const
         }
     }
     for (size_t k = 0; k < m_kk; k++) {
-        m_tmpV[k] = 0.0;
+        m_workS[k] = 0.0;
         for (size_t i = 0; i < m_kk; i++) {
             size_t counter = k + m_kk*i;
-            m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
+            m_workS[k] += moleFractions_[i] * a_coeff_vec(1,counter);
         }
     }
 
@@ -293,7 +293,7 @@ void RedlichKwongMFTP::getPartialMolarEntropies(double* sbar) const
                    + GasConstant * log(mv/vmb)
                    + GasConstant * b_vec_Curr_[k]/vmb
                    + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
-                   - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
+                   - 2.0 * m_workS[k]/(m_b_current * sqt) * log(vpb/mv)
                    + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
                    - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
                   );
@@ -327,10 +327,10 @@ void RedlichKwongMFTP::getPartialMolarVolumes(double* vbar) const
         }
     }
     for (size_t k = 0; k < m_kk; k++) {
-        m_tmpV[k] = 0.0;
+        m_workS[k] = 0.0;
         for (size_t i = 0; i < m_kk; i++) {
             size_t counter = k + m_kk*i;
-            m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
+            m_workS[k] += moleFractions_[i] * a_coeff_vec(1,counter);
         }
     }
 
diff --git a/src/thermo/VPStandardStateTP.cpp b/src/thermo/VPStandardStateTP.cpp
index 26dd080ef2..37d7e1914b 100644
--- a/src/thermo/VPStandardStateTP.cpp
+++ b/src/thermo/VPStandardStateTP.cpp
@@ -199,8 +199,8 @@ void VPStandardStateTP::setPressure(double p)
 
 void VPStandardStateTP::calcDensity()
 {
-    getPartialMolarVolumes(m_tmpV.data());
-    double dd = meanMolecularWeight() / mean_X(m_tmpV);
+    getPartialMolarVolumes(m_workS.data());
+    double dd = meanMolecularWeight() / mean_X(m_workS);
     Phase::assignDensity(dd);
 }