diff --git a/include/cantera/thermo/ConstDensityThermo.h b/include/cantera/thermo/ConstDensityThermo.h
index d89353bcb..5153e572a 100644
--- a/include/cantera/thermo/ConstDensityThermo.h
+++ b/include/cantera/thermo/ConstDensityThermo.h
@@ -121,7 +121,7 @@ public:
         return m_cp0_R;
     }
 
-    virtual void initThermo();
+    virtual bool addSpecies(shared_ptr<Species> spec);
 
     virtual void setParameters(int n, doublereal* const c) {
         setDensity(c[0]);
diff --git a/include/cantera/thermo/DebyeHuckel.h b/include/cantera/thermo/DebyeHuckel.h
index 869e3fca2..18e21aff6 100644
--- a/include/cantera/thermo/DebyeHuckel.h
+++ b/include/cantera/thermo/DebyeHuckel.h
@@ -829,7 +829,7 @@ public:
      *  -------------- Utilities -------------------------------
      */
 
-    virtual void initThermo();
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
 
     //! Return the Debye Huckel constant as a function of temperature
@@ -1112,14 +1112,6 @@ protected:
      */
     vector_fp m_B_Dot;
 
-    /**
-     *  These are coefficients to describe the increase in activity coeff for
-     *  non-polar molecules due to the electrolyte becoming stronger (the
-     *  so-called salt-out effect)
-     */
-    vector_fp m_npActCoeff;
-
-
     //! Pointer to the Water standard state object
     /*!
      *  derived from the equation of state for water.
@@ -1180,13 +1172,6 @@ protected:
     mutable vector_fp m_dlnActCoeffMolaldP;
 
 private:
-    //! Initialize the internal lengths.
-    /*!
-     * This internal function adjusts the lengths of arrays based on
-     * the number of species.
-     */
-    void initLengths();
-
     //! Calculate the log activity coefficients
     /*!
      * This function updates the internally stored natural logarithm of the
diff --git a/include/cantera/thermo/GibbsExcessVPSSTP.h b/include/cantera/thermo/GibbsExcessVPSSTP.h
index 4bf61d420..4520275ac 100644
--- a/include/cantera/thermo/GibbsExcessVPSSTP.h
+++ b/include/cantera/thermo/GibbsExcessVPSSTP.h
@@ -244,12 +244,7 @@ public:
     virtual void setConcentrations(const doublereal* const c);
     //@}
 
-    virtual void initThermo();
-
-private:
-    //! Initialize lengths of local variables after all species have
-    //! been identified.
-    void initLengths();
+    virtual bool addSpecies(shared_ptr<Species> spec);
 
 protected:
     //! utility routine to check mole fraction sum
diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h
index 23ac2e799..f3f8c1e04 100644
--- a/include/cantera/thermo/IdealGasPhase.h
+++ b/include/cantera/thermo/IdealGasPhase.h
@@ -605,7 +605,7 @@ public:
 
     //@}
 
-    virtual void initThermo();
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void setToEquilState(const doublereal* lambda_RT);
 
 protected:
diff --git a/include/cantera/thermo/IdealMolalSoln.h b/include/cantera/thermo/IdealMolalSoln.h
index 3e6906da9..b5f81ebd7 100644
--- a/include/cantera/thermo/IdealMolalSoln.h
+++ b/include/cantera/thermo/IdealMolalSoln.h
@@ -412,7 +412,7 @@ public:
 
     // -------------- Utilities -------------------------------
 
-    virtual void initThermo();
+    virtual bool addSpecies(shared_ptr<Species> spec);
 
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id="");
 
@@ -513,12 +513,6 @@ private:
      */
     void s_updateIMS_lnMolalityActCoeff() const;
 
-    //! This internal function adjusts the lengths of arrays.
-    /*!
-     * This function is not virtual nor is it inherited
-     */
-    void initLengths();
-
     //! Calculate parameters for cutoff treatments of activity coefficients
     /*!
      * Some cutoff treatments for the activity coefficients actually require
diff --git a/include/cantera/thermo/IdealSolidSolnPhase.h b/include/cantera/thermo/IdealSolidSolnPhase.h
index 8ef4eb28f..eaf269784 100644
--- a/include/cantera/thermo/IdealSolidSolnPhase.h
+++ b/include/cantera/thermo/IdealSolidSolnPhase.h
@@ -619,8 +619,8 @@ public:
     /// @name Utility Functions
     //@{
 
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
-
     virtual void setToEquilState(const doublereal* lambda_RT);
 
     /**
@@ -718,9 +718,6 @@ private:
      */
     void _updateThermo() const;
 
-    //! This internal function adjusts the lengths of arrays
-    void initLengths();
-
     //@}
 };
 }
diff --git a/include/cantera/thermo/IdealSolnGasVPSS.h b/include/cantera/thermo/IdealSolnGasVPSS.h
index bf0203d0d..c097e6525 100644
--- a/include/cantera/thermo/IdealSolnGasVPSS.h
+++ b/include/cantera/thermo/IdealSolnGasVPSS.h
@@ -144,18 +144,11 @@ public:
      */
     //@{
 
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void setParametersFromXML(const XML_Node& thermoNode);
-    virtual void initThermo();
     virtual void setToEquilState(const doublereal* lambda_RT);
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
 
-private:
-    //!  @internal Initialize the internal lengths in this object.
-    /*!
-     * Note this is not a virtual function and only handles this object
-     */
-    void initLengths();
-
     //@}
 
 protected:
diff --git a/include/cantera/thermo/LatticePhase.h b/include/cantera/thermo/LatticePhase.h
index bcd4f2928..e7f165cad 100644
--- a/include/cantera/thermo/LatticePhase.h
+++ b/include/cantera/thermo/LatticePhase.h
@@ -625,7 +625,7 @@ public:
     /// @name  Utilities for Initialization of the Object
     //@{
 
-    virtual void initThermo();
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
 
     //! Set the equation of state parameters from the argument list
diff --git a/include/cantera/thermo/LatticeSolidPhase.h b/include/cantera/thermo/LatticeSolidPhase.h
index f0bf0d56e..af237589b 100644
--- a/include/cantera/thermo/LatticeSolidPhase.h
+++ b/include/cantera/thermo/LatticeSolidPhase.h
@@ -421,11 +421,9 @@ public:
     virtual void getGibbs_RT_ref(doublereal* grt) const;
     virtual void getGibbs_ref(doublereal* g) const;
 
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void initThermo();
 
-    //! Initialize vectors that depend on the number of species and sublattices
-    void initLengths();
-
     //! Add in species from Slave phases
     /*!
      *  This hook is used for cSS_CONVENTION_SLAVE phases
diff --git a/include/cantera/thermo/MixtureFugacityTP.h b/include/cantera/thermo/MixtureFugacityTP.h
index 2b4932f79..bdea15983 100644
--- a/include/cantera/thermo/MixtureFugacityTP.h
+++ b/include/cantera/thermo/MixtureFugacityTP.h
@@ -372,14 +372,8 @@ public:
      */
     //@{
 
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void setStateFromXML(const XML_Node& state);
-    virtual void initThermo();
-    virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
-
-private:
-    //!  @internal Initialize the internal lengths in this object.
-    void initLengths();
-    //@}
 
 protected:
     //! @name Special Functions for fugacity classes
diff --git a/include/cantera/thermo/MolalityVPSSTP.h b/include/cantera/thermo/MolalityVPSSTP.h
index 8a914fb27..6ad2d99c0 100644
--- a/include/cantera/thermo/MolalityVPSSTP.h
+++ b/include/cantera/thermo/MolalityVPSSTP.h
@@ -485,8 +485,8 @@ public:
     /// see importPhase().
     //@{
 
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void initThermo();
-    virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
 
     //@}
 
@@ -567,10 +567,6 @@ private:
      */
     virtual size_t findCLMIndex() const;
 
-    //! Initialize lengths of local variables after all species have been
-    //! identified.
-    void initLengths();
-
 protected:
 
     //! Index of the solvent. Currently the index of the solvent is hard-coded
diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h
index 96eb1827d..180cb8d13 100644
--- a/include/cantera/thermo/Phase.h
+++ b/include/cantera/thermo/Phase.h
@@ -686,6 +686,13 @@ public:
 
     //! Add a Species to this Phase. Returns `true` if the species was
     //! successfully added, or `false` if the species was ignored.
+    //!
+    //! Derived classes which need to size arrays according to the number of
+    //! species should overload this method. The derived class implementation
+    //! should call the base class method, and, if this returns `true`
+    //! (indicating that the species has been added), adjust their array sizes
+    //! accordingly.
+    //!
     //! @see ignoreUndefinedElements addUndefinedElements throwUndefinedElements
     virtual bool addSpecies(shared_ptr<Species> spec);
 
diff --git a/include/cantera/thermo/RedlichKwongMFTP.h b/include/cantera/thermo/RedlichKwongMFTP.h
index 688660452..b7d3d09e3 100644
--- a/include/cantera/thermo/RedlichKwongMFTP.h
+++ b/include/cantera/thermo/RedlichKwongMFTP.h
@@ -190,8 +190,8 @@ public:
      */
     //@{
 
+    virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void setParametersFromXML(const XML_Node& thermoNode);
-    virtual void initThermo();
     virtual void setToEquilState(const doublereal* lambda_RT);
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
 
@@ -211,11 +211,6 @@ private:
      */
     void readXMLCrossFluid(XML_Node& pureFluidParam);
 
-    //!  @internal Initialize the internal lengths in this object.
-    /*!
-     * Note this is not a virtual function and only handles this object
-     */
-    void initLengths();
     // @}
 
 protected:
diff --git a/include/cantera/thermo/SingleSpeciesTP.h b/include/cantera/thermo/SingleSpeciesTP.h
index 32ed6dd0f..c5519461e 100644
--- a/include/cantera/thermo/SingleSpeciesTP.h
+++ b/include/cantera/thermo/SingleSpeciesTP.h
@@ -252,13 +252,7 @@ public:
                              doublereal tol = 1.e-8);
     //@}
 
-    /**
-     * @copydoc ThermoPhase::initThermo
-     *
-     * This version sets the mole fraction vector to x[0] = 1.0, and then
-     * calls the ThermoPhase::initThermo() function.
-     */
-    virtual void initThermo();
+    virtual bool addSpecies(shared_ptr<Species> spec);
 
 protected:
     //! The current pressure of the solution (Pa). It gets initialized to 1 atm.
diff --git a/include/cantera/thermo/SurfPhase.h b/include/cantera/thermo/SurfPhase.h
index b5053b2dc..8c1c14ba0 100644
--- a/include/cantera/thermo/SurfPhase.h
+++ b/include/cantera/thermo/SurfPhase.h
@@ -293,8 +293,7 @@ public:
      */
     virtual void setParametersFromXML(const XML_Node& thermoData);
 
-    //! Initialize the SurfPhase object after all species have been set up
-    virtual void initThermo();
+    virtual bool addSpecies(shared_ptr<Species> spec);
 
     //! Set the initial state of the Surface Phase from an XML_Node
     /*!
diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h
index ab726e3c3..e3f4eec9c 100644
--- a/include/cantera/thermo/ThermoPhase.h
+++ b/include/cantera/thermo/ThermoPhase.h
@@ -1381,7 +1381,6 @@ public:
      */
     //@{
 
-    using Phase::addSpecies;
     virtual bool addSpecies(shared_ptr<Species> spec);
 
     virtual void modifySpecies(size_t k, shared_ptr<Species> spec);
diff --git a/src/thermo/ConstDensityThermo.cpp b/src/thermo/ConstDensityThermo.cpp
index 96e9abff2..96d8a2aec 100644
--- a/src/thermo/ConstDensityThermo.cpp
+++ b/src/thermo/ConstDensityThermo.cpp
@@ -109,14 +109,17 @@ void ConstDensityThermo::getStandardChemPotentials(doublereal* mu0) const
     getPureGibbs(mu0);
 }
 
-void ConstDensityThermo::initThermo()
+bool ConstDensityThermo::addSpecies(shared_ptr<Species> spec)
 {
-    ThermoPhase::initThermo();
-    m_h0_RT.resize(m_kk);
-    m_g0_RT.resize(m_kk);
-    m_cp0_R.resize(m_kk);
-    m_s0_R.resize(m_kk);
-    m_pp.resize(m_kk);
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        m_h0_RT.push_back(0.0);
+        m_g0_RT.push_back(0.0);
+        m_cp0_R.push_back(0.0);
+        m_s0_R.push_back(0.0);
+        m_pp.push_back(0.0);
+    }
+    return added;
 }
 
 void ConstDensityThermo::_updateThermo() const
diff --git a/src/thermo/DebyeHuckel.cpp b/src/thermo/DebyeHuckel.cpp
index 7c283b52c..2e9eada3b 100644
--- a/src/thermo/DebyeHuckel.cpp
+++ b/src/thermo/DebyeHuckel.cpp
@@ -100,7 +100,6 @@ DebyeHuckel& DebyeHuckel::operator=(const DebyeHuckel& b)
         m_formDH = b.m_formDH;
         m_formGC = b.m_formGC;
         m_Aionic = b.m_Aionic;
-        m_npActCoeff = b.m_npActCoeff;
         m_IionicMolality = b.m_IionicMolality;
         m_maxIionicStrength = b.m_maxIionicStrength;
         m_useHelgesonFixedForm= b.m_useHelgesonFixedForm;
@@ -109,7 +108,6 @@ DebyeHuckel& DebyeHuckel::operator=(const DebyeHuckel& b)
         m_A_Debye = b.m_A_Debye;
         m_B_Debye = b.m_B_Debye;
         m_B_Dot = b.m_B_Dot;
-        m_npActCoeff = b.m_npActCoeff;
 
         // This is an internal shallow copy of the PDSS_Water pointer
         m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
@@ -400,16 +398,6 @@ void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const
 
 // -------------- Utilities -------------------------------
 
-void DebyeHuckel::initThermo()
-{
-    MolalityVPSSTP::initThermo();
-    m_npActCoeff.resize(3);
-    m_npActCoeff[0] = 0.1127;
-    m_npActCoeff[1] = -0.01049;
-    m_npActCoeff[2] = 1.545E-3;
-    initLengths();
-}
-
 //! Utility function to assign an integer value from a string for the
 //! ElectrolyteSpeciesType field.
 /*!
@@ -542,10 +530,6 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
                            " should be first species");
     }
 
-    // Initialize all of the lengths of arrays in the object now that we know
-    // what species are in the phase.
-    initThermo();
-
     // Now go get the specification of the standard states for species in the
     // solution. This includes the molar volumes data blocks for incompressible
     // species.
@@ -727,6 +711,7 @@ void DebyeHuckel::initThermoXML(XML_Node& phaseNode, const std::string& id_)
         if (acNode.hasChild("DHBetaMatrix")) {
             if (m_formDH == DHFORM_BETAIJ ||
                     m_formDH == DHFORM_PITZER_BETAIJ) {
+                m_Beta_ij.resize(m_kk, m_kk, 0.0);
                 XML_Node& irNode = acNode.child("DHBetaMatrix");
                 const vector<string>& sn = speciesNames();
                 getMatrixValues(irNode, sn, sn, m_Beta_ij, true, true);
@@ -942,31 +927,35 @@ double DebyeHuckel::AionicRadius(int k) const
 
 // ------------ Private and Restricted Functions ------------------
 
-void DebyeHuckel::initLengths()
+bool DebyeHuckel::addSpecies(shared_ptr<Species> spec)
 {
-    m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
-    m_speciesSize.resize(m_kk);
-    m_Aionic.resize(m_kk, 0.0);
-    m_lnActCoeffMolal.resize(m_kk, 0.0);
-    m_dlnActCoeffMolaldT.resize(m_kk, 0.0);
-    m_d2lnActCoeffMolaldT2.resize(m_kk, 0.0);
-    m_dlnActCoeffMolaldP.resize(m_kk, 0.0);
-    m_B_Dot.resize(m_kk, 0.0);
-    m_pp.resize(m_kk, 0.0);
-    m_tmpV.resize(m_kk, 0.0);
-    if (m_formDH == DHFORM_BETAIJ ||
-            m_formDH == DHFORM_PITZER_BETAIJ) {
-        m_Beta_ij.resize(m_kk, m_kk, 0.0);
+    bool added = MolalityVPSSTP::addSpecies(spec);
+    if (added) {
+        m_electrolyteSpeciesType.push_back(cEST_polarNeutral);
+        m_speciesSize.push_back(0.0);
+        m_Aionic.push_back(0.0);
+        m_lnActCoeffMolal.push_back(0.0);
+        m_dlnActCoeffMolaldT.push_back(0.0);
+        m_d2lnActCoeffMolaldT2.push_back(0.0);
+        m_dlnActCoeffMolaldP.push_back(0.0);
+        m_B_Dot.push_back(0.0);
+        m_pp.push_back(0.0);
+        m_tmpV.push_back(0.0);
     }
+    return added;
 }
 
 double DebyeHuckel::_nonpolarActCoeff(double IionicMolality) const
 {
+     // These are coefficients to describe the increase in activity coeff for
+     // non-polar molecules due to the electrolyte becoming stronger (the so-
+     // called salt-out effect)
+    const static double npActCoeff[] = {0.1127, -0.01049, 1.545E-3};
     double I2 = IionicMolality * IionicMolality;
     double l10actCoeff =
-        m_npActCoeff[0] * IionicMolality +
-        m_npActCoeff[1] * I2 +
-        m_npActCoeff[2] * I2 * IionicMolality;
+        npActCoeff[0] * IionicMolality +
+        npActCoeff[1] * I2 +
+        npActCoeff[2] * I2 * IionicMolality;
     return pow(10.0 , l10actCoeff);
 }
 
diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp
index 1543a612d..683a17c03 100644
--- a/src/thermo/GibbsExcessVPSSTP.cpp
+++ b/src/thermo/GibbsExcessVPSSTP.cpp
@@ -157,23 +157,24 @@ double GibbsExcessVPSSTP::checkMFSum(const doublereal* const x) const
     return norm;
 }
 
-void GibbsExcessVPSSTP::initThermo()
+bool GibbsExcessVPSSTP::addSpecies(shared_ptr<Species> spec)
 {
-    initLengths();
-    VPStandardStateTP::initThermo();
-    getMoleFractions(moleFractions_.data());
-}
-
-void GibbsExcessVPSSTP::initLengths()
-{
-    moleFractions_.resize(m_kk);
-    lnActCoeff_Scaled_.resize(m_kk);
-    dlnActCoeffdT_Scaled_.resize(m_kk);
-    d2lnActCoeffdT2_Scaled_.resize(m_kk);
-    dlnActCoeffdlnX_diag_.resize(m_kk);
-    dlnActCoeffdlnN_diag_.resize(m_kk);
-    dlnActCoeffdlnN_.resize(m_kk, m_kk);
-    m_pp.resize(m_kk);
+    bool added = VPStandardStateTP::addSpecies(spec);
+    if (added) {
+        if (m_kk == 1) {
+            moleFractions_.push_back(1.0);
+        } else {
+            moleFractions_.push_back(0.0);
+        }
+        lnActCoeff_Scaled_.push_back(0.0);
+        dlnActCoeffdT_Scaled_.push_back(0.0);
+        d2lnActCoeffdT2_Scaled_.push_back(0.0);
+        dlnActCoeffdlnX_diag_.push_back(0.0);
+        dlnActCoeffdlnN_diag_.push_back(0.0);
+        dlnActCoeffdlnN_.resize(m_kk, m_kk);
+        m_pp.push_back(0.0);
+    }
+    return added;
 }
 
 } // end of namespace Cantera
diff --git a/src/thermo/IdealGasPhase.cpp b/src/thermo/IdealGasPhase.cpp
index ce0ac3990..96a0a48ca 100644
--- a/src/thermo/IdealGasPhase.cpp
+++ b/src/thermo/IdealGasPhase.cpp
@@ -261,16 +261,21 @@ void IdealGasPhase::getStandardVolumes_ref(doublereal* vol) const
     }
 }
 
-void IdealGasPhase::initThermo()
+bool IdealGasPhase::addSpecies(shared_ptr<Species> spec)
 {
-    ThermoPhase::initThermo();
-    m_p0 = refPressure();
-    m_h0_RT.resize(m_kk);
-    m_g0_RT.resize(m_kk);
-    m_expg0_RT.resize(m_kk);
-    m_cp0_R.resize(m_kk);
-    m_s0_R.resize(m_kk);
-    m_pp.resize(m_kk);
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        if (m_kk == 1) {
+            m_p0 = refPressure();
+        }
+        m_h0_RT.push_back(0.0);
+        m_g0_RT.push_back(0.0);
+        m_expg0_RT.push_back(0.0);
+        m_cp0_R.push_back(0.0);
+        m_s0_R.push_back(0.0);
+        m_pp.push_back(0.0);
+    }
+    return added;
 }
 
 void IdealGasPhase::setToEquilState(const doublereal* mu_RT)
diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp
index 2a054f187..60e392c77 100644
--- a/src/thermo/IdealMolalSoln.cpp
+++ b/src/thermo/IdealMolalSoln.cpp
@@ -393,10 +393,16 @@ void IdealMolalSoln::getPartialMolarCp(doublereal* cpbar) const
 
 // -------------- Utilities -------------------------------
 
-void IdealMolalSoln::initThermo()
+bool IdealMolalSoln::addSpecies(shared_ptr<Species> spec)
 {
-    initLengths();
-    MolalityVPSSTP::initThermo();
+    bool added = MolalityVPSSTP::addSpecies(spec);
+    if (added) {
+        m_pp.push_back(0.0);
+        m_speciesMolarVolume.push_back(0.0);
+        m_tmpV.push_back(0.0);
+        IMS_lnActCoeffMolal_.push_back(0.0);
+    }
+    return added;
 }
 
 void IdealMolalSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
@@ -641,14 +647,6 @@ void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
     }
 }
 
-void IdealMolalSoln::initLengths()
-{
-    m_pp.resize(m_kk);
-    m_speciesMolarVolume.resize(m_kk);
-    m_tmpV.resize(m_kk);
-    IMS_lnActCoeffMolal_.resize(m_kk);
-}
-
 void IdealMolalSoln::calcIMSCutoffParams_()
 {
     IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
diff --git a/src/thermo/IdealSolidSolnPhase.cpp b/src/thermo/IdealSolidSolnPhase.cpp
index 685a05718..060dfb485 100644
--- a/src/thermo/IdealSolidSolnPhase.cpp
+++ b/src/thermo/IdealSolidSolnPhase.cpp
@@ -450,6 +450,30 @@ const vector_fp& IdealSolidSolnPhase::entropy_R_ref() const
 
 // Utility Functions
 
+bool IdealSolidSolnPhase::addSpecies(shared_ptr<Species> spec)
+{
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        if (m_kk == 1) {
+            // Obtain the reference pressure by calling the ThermoPhase function
+            // refPressure, which in turn calls the species thermo reference
+            // pressure function of the same name.
+            m_Pref = refPressure();
+        }
+
+        m_h0_RT.push_back(0.0);
+        m_g0_RT.push_back(0.0);
+        m_expg0_RT.push_back(0.0);
+        m_cp0_R.push_back(0.0);
+        m_s0_R.push_back(0.0);
+        m_pe.push_back(0.0);;
+        m_pp.push_back(0.0);
+        m_speciesMolarVolume.push_back(0.0);
+    }
+    return added;
+}
+
+
 void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
 {
     if (id_.size() > 0 && phaseNode.id() != id_) {
@@ -495,10 +519,6 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
                            "Unspecified standardConc model");
     }
 
-    // Initialize all of the lengths now that we know how many species
-    // there are in the phase.
-    initLengths();
-
     // Now go get the molar volumes
     XML_Node& speciesList = phaseNode.child("speciesArray");
     XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"],
@@ -514,23 +534,6 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
     ThermoPhase::initThermoXML(phaseNode, id_);
 }
 
-void IdealSolidSolnPhase::initLengths()
-{
-    // Obtain the reference pressure by calling the ThermoPhase function
-    // refPressure, which in turn calls the species thermo reference pressure
-    // function of the same name.
-    m_Pref = refPressure();
-
-    m_h0_RT.resize(m_kk);
-    m_g0_RT.resize(m_kk);
-    m_expg0_RT.resize(m_kk);
-    m_cp0_R.resize(m_kk);
-    m_s0_R.resize(m_kk);
-    m_pe.resize(m_kk, 0.0);
-    m_pp.resize(m_kk);
-    m_speciesMolarVolume.resize(m_kk);
-}
-
 void IdealSolidSolnPhase::setToEquilState(const doublereal* lambda_RT)
 {
     const vector_fp& grt = gibbs_RT_ref();
diff --git a/src/thermo/IdealSolnGasVPSS.cpp b/src/thermo/IdealSolnGasVPSS.cpp
index ee636f28a..f85245687 100644
--- a/src/thermo/IdealSolnGasVPSS.cpp
+++ b/src/thermo/IdealSolnGasVPSS.cpp
@@ -233,12 +233,6 @@ void IdealSolnGasVPSS::getPartialMolarVolumes(doublereal* vbar) const
     getStandardVolumes(vbar);
 }
 
-void IdealSolnGasVPSS::initThermo()
-{
-    initLengths();
-    VPStandardStateTP::initThermo();
-}
-
 void IdealSolnGasVPSS::setToEquilState(const doublereal* mu_RT)
 {
     updateStandardStateThermo();
@@ -268,15 +262,17 @@ void IdealSolnGasVPSS::setToEquilState(const doublereal* mu_RT)
     setState_PX(pres, &m_pp[0]);
 }
 
-void IdealSolnGasVPSS::initLengths()
+bool IdealSolnGasVPSS::addSpecies(shared_ptr<Species> spec)
 {
-    m_pp.resize(m_kk, 0.0);
+    bool added = VPStandardStateTP::addSpecies(spec);
+    if (added) {
+        m_pp.push_back(0.0);
+    }
+    return added;
 }
 
 void IdealSolnGasVPSS::initThermoXML(XML_Node& phaseNode, const std::string& id_)
 {
-    IdealSolnGasVPSS::initLengths();
-
     if (phaseNode.hasChild("thermo")) {
         XML_Node& thermoNode = phaseNode.child("thermo");
         std::string model = thermoNode["model"];
diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp
index e3e26889d..dd1974a22 100644
--- a/src/thermo/LatticePhase.cpp
+++ b/src/thermo/LatticePhase.cpp
@@ -277,16 +277,20 @@ const vector_fp& LatticePhase::cp_R_ref() const
     return m_cp0_R;
 }
 
-void LatticePhase::initThermo()
+bool LatticePhase::addSpecies(shared_ptr<Species> spec)
 {
-    m_Pref = refPressure();
-    m_h0_RT.resize(m_kk);
-    m_g0_RT.resize(m_kk);
-    m_cp0_R.resize(m_kk);
-    m_s0_R.resize(m_kk);
-    m_speciesMolarVolume.resize(m_kk, 0.0);
-
-    ThermoPhase::initThermo();
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        if (m_kk == 1) {
+            m_Pref = refPressure();
+        }
+        m_h0_RT.push_back(0.0);
+        m_g0_RT.push_back(0.0);
+        m_cp0_R.push_back(0.0);
+        m_s0_R.push_back(0.0);
+        m_speciesMolarVolume.push_back(0.0);
+    }
+    return added;
 }
 
 void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id_)
diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp
index 3b8080350..d91e86709 100644
--- a/src/thermo/LatticeSolidPhase.cpp
+++ b/src/thermo/LatticeSolidPhase.cpp
@@ -360,7 +360,7 @@ void LatticeSolidPhase::installSlavePhases(XML_Node* phaseNode)
 
 void LatticeSolidPhase::initThermo()
 {
-    initLengths();
+    lkstart_.resize(m_lattice.size() + 1);
     size_t loc = 0;
     for (size_t n = 0; n < m_lattice.size(); n++) {
         size_t nsp = m_lattice[n]->nSpecies();
@@ -375,12 +375,14 @@ void LatticeSolidPhase::initThermo()
     ThermoPhase::initThermo();
 }
 
-void LatticeSolidPhase::initLengths()
+bool LatticeSolidPhase::addSpecies(shared_ptr<Species> spec)
 {
-    theta_.resize(m_lattice.size(), 0);
-    lkstart_.resize(m_lattice.size() + 1);
-    m_x.resize(m_kk, 0.0);
-    tmpV_.resize(m_kk, 0.0);
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        m_x.push_back(0.0);
+        tmpV_.push_back(0.0);
+    }
+    return added;
 }
 
 void LatticeSolidPhase::_updateThermo() const
diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp
index e1110001d..f1510a57d 100644
--- a/src/thermo/MixtureFugacityTP.cpp
+++ b/src/thermo/MixtureFugacityTP.cpp
@@ -249,20 +249,21 @@ void MixtureFugacityTP::setStateFromXML(const XML_Node& state)
     }
 }
 
-void MixtureFugacityTP::initThermo()
+bool MixtureFugacityTP::addSpecies(shared_ptr<Species> spec)
 {
-    initLengths();
-    ThermoPhase::initThermo();
-}
-
-void MixtureFugacityTP::initLengths()
-{
-    moleFractions_.resize(m_kk, 0.0);
-    moleFractions_[0] = 1.0;
-    m_h0_RT.resize(m_kk, 0.0);
-    m_cp0_R.resize(m_kk, 0.0);
-    m_g0_RT.resize(m_kk, 0.0);
-    m_s0_R.resize(m_kk, 0.0);
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        if (m_kk == 1) {
+            moleFractions_.push_back(1.0);
+        } else {
+            moleFractions_.push_back(0.0);
+        }
+        m_h0_RT.push_back(0.0);
+        m_cp0_R.push_back(0.0);
+        m_g0_RT.push_back(0.0);
+        m_s0_R.push_back(0.0);
+    }
+    return added;
 }
 
 void MixtureFugacityTP::setTemperature(const doublereal temp)
@@ -409,12 +410,6 @@ void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublerea
     setState_TP(t,p);
 }
 
-void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
-{
-    MixtureFugacityTP::initLengths();
-    ThermoPhase::initThermoXML(phaseNode, id_);
-}
-
 doublereal MixtureFugacityTP::z() const
 {
     return pressure() * meanMolecularWeight() / (density() * RT());
diff --git a/src/thermo/MolalityVPSSTP.cpp b/src/thermo/MolalityVPSSTP.cpp
index 166fdb0c6..cbadc75d8 100644
--- a/src/thermo/MolalityVPSSTP.cpp
+++ b/src/thermo/MolalityVPSSTP.cpp
@@ -325,12 +325,8 @@ void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, const std::string&
 
 void MolalityVPSSTP::initThermo()
 {
-    initLengths();
     VPStandardStateTP::initThermo();
 
-    // The solvent defaults to species 0
-    setSolvent(0);
-
     // Find the Cl- species
     m_indexCLM = findCLMIndex();
 }
@@ -401,20 +397,17 @@ size_t MolalityVPSSTP::findCLMIndex() const
     return indexCLM;
 }
 
-//   Initialize lengths of local variables after all species have
-//   been identified.
-void MolalityVPSSTP::initLengths()
+bool MolalityVPSSTP::addSpecies(shared_ptr<Species> spec)
 {
-    m_molalities.resize(m_kk);
-}
-
-void MolalityVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
-{
-    initLengths();
-
-    // The solvent defaults to species 0
-    setSolvent(0);
-    VPStandardStateTP::initThermoXML(phaseNode, id_);
+    bool added = VPStandardStateTP::addSpecies(spec);
+    if (added) {
+        if (m_kk == 1) {
+            // The solvent defaults to species 0
+            setSolvent(0);
+        }
+        m_molalities.push_back(0.0);
+    }
+    return added;
 }
 
 std::string MolalityVPSSTP::report(bool show_thermo, doublereal threshold) const
diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp
index ded087320..de071e5d4 100644
--- a/src/thermo/Phase.cpp
+++ b/src/thermo/Phase.cpp
@@ -834,6 +834,7 @@ bool Phase::addSpecies(shared_ptr<Species> spec) {
         m_y.push_back(0.0);
         m_ym.push_back(0.0);
     }
+    invalidateCache();
     return true;
 }
 
diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp
index 56b03bec3..3c3ee3107 100644
--- a/src/thermo/RedlichKwongMFTP.cpp
+++ b/src/thermo/RedlichKwongMFTP.cpp
@@ -34,6 +34,7 @@ RedlichKwongMFTP::RedlichKwongMFTP() :
     m_b_current(0.0),
     m_a_current(0.0),
     NSolns_(0),
+    Vroot_{0.0, 0.0, 0.0},
     dpdV_(0.0),
     dpdT_(0.0)
 {
@@ -45,6 +46,7 @@ RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, const std::string&
     m_b_current(0.0),
     m_a_current(0.0),
     NSolns_(0),
+    Vroot_{0.0, 0.0, 0.0},
     dpdV_(0.0),
     dpdT_(0.0)
 {
@@ -57,6 +59,7 @@ RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id
     m_b_current(0.0),
     m_a_current(0.0),
     NSolns_(0),
+    Vroot_{0.0, 0.0, 0.0},
     dpdV_(0.0),
     dpdT_(0.0)
 {
@@ -541,16 +544,6 @@ doublereal RedlichKwongMFTP::critDensity() const
     return mmw / vc;
 }
 
-void RedlichKwongMFTP::initThermo()
-{
-    Vroot_[0] = 0.0;
-    Vroot_[1] = 0.0;
-    Vroot_[2] = 0.0;
-
-    initLengths();
-    MixtureFugacityTP::initThermo();
-}
-
 void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
 {
     double tmp, tmp2;
@@ -581,27 +574,29 @@ void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
     setState_PX(pres, &m_pp[0]);
 }
 
-void RedlichKwongMFTP::initLengths()
+bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
 {
-    a_vec_Curr_.resize(m_kk * m_kk, 0.0);
-    b_vec_Curr_.resize(m_kk, 0.0);
+    bool added = MixtureFugacityTP::addSpecies(spec);
+    if (added) {
+        a_vec_Curr_.resize(m_kk * m_kk, 0.0);
+        b_vec_Curr_.push_back(0.0);
 
-    a_coeff_vec.resize(2, m_kk * m_kk, 0.0);
+        a_coeff_vec.resize(2, m_kk * m_kk, 0.0);
 
-    m_pc_Species.resize(m_kk, 0.0);
-    m_tc_Species.resize(m_kk, 0.0);
-    m_vc_Species.resize(m_kk, 0.0);
+        m_pc_Species.push_back(0.0);
+        m_tc_Species.push_back(0.0);
+        m_vc_Species.push_back(0.0);
 
-    m_pp.resize(m_kk, 0.0);
-    m_tmpV.resize(m_kk, 0.0);
-    m_partialMolarVolumes.resize(m_kk, 0.0);
-    dpdni_.resize(m_kk, 0.0);
+        m_pp.push_back(0.0);
+        m_tmpV.push_back(0.0);
+        m_partialMolarVolumes.push_back(0.0);
+        dpdni_.push_back(0.0);
+    }
+    return added;
 }
 
 void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
 {
-    RedlichKwongMFTP::initLengths();
-
     // Check the model parameter for the Redlich-Kwong equation of state
     // two are allowed
     //       RedlichKwong        mixture of species, each of which are RK fluids
diff --git a/src/thermo/SingleSpeciesTP.cpp b/src/thermo/SingleSpeciesTP.cpp
index ccb7ebfe7..a57f476bb 100644
--- a/src/thermo/SingleSpeciesTP.cpp
+++ b/src/thermo/SingleSpeciesTP.cpp
@@ -20,7 +20,10 @@ namespace Cantera
 {
 SingleSpeciesTP::SingleSpeciesTP() :
     m_press(OneAtm),
-    m_p0(OneAtm)
+    m_p0(OneAtm),
+    m_h0_RT(1),
+    m_cp0_R(1),
+    m_s0_R(1)
 {
 }
 
@@ -276,25 +279,18 @@ void SingleSpeciesTP::setState_SV(doublereal s, doublereal v,
     throw CanteraError("setState_SV","no convergence. dt = {}", dt);
 }
 
-void SingleSpeciesTP::initThermo()
+bool SingleSpeciesTP::addSpecies(shared_ptr<Species> spec)
 {
-    // Make sure there is one and only one species in this phase.
-    if (nSpecies() != 1) {
-        throw CanteraError("initThermo",
-                           "stoichiometric substances may only contain one species.");
+    if (m_kk != 0) {
+        throw CanteraError("SingleSpeciesTP::addSpecies",
+            "Stoichiometric substances may only contain one species.");
     }
-
-    // Resize temporary arrays.
-    m_h0_RT.resize(1);
-    m_cp0_R.resize(1);
-    m_s0_R.resize(1);
-
-    // Make sure the species mole fraction is equal to 1.0;
-    double x = 1.0;
-    ThermoPhase::setMoleFractions(&x);
-
-    // Call the base class initThermo object.
-    ThermoPhase::initThermo();
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        double x = 1.0;
+        ThermoPhase::setMoleFractions(&x);
+    }
+    return added;
 }
 
 void SingleSpeciesTP::_updateThermo() const
diff --git a/src/thermo/SurfPhase.cpp b/src/thermo/SurfPhase.cpp
index c5f9f8794..bcdfab9e2 100644
--- a/src/thermo/SurfPhase.cpp
+++ b/src/thermo/SurfPhase.cpp
@@ -232,21 +232,22 @@ void SurfPhase::getCp_R_ref(doublereal* cprt) const
     getCp_R(cprt);
 }
 
-void SurfPhase::initThermo()
+bool SurfPhase::addSpecies(shared_ptr<Species> spec)
 {
-    ThermoPhase::initThermo();
-    m_h0.resize(m_kk);
-    m_s0.resize(m_kk);
-    m_cp0.resize(m_kk);
-    m_mu0.resize(m_kk);
-    m_work.resize(m_kk);
-    vector_fp cov(m_kk, 0.0);
-    cov[0] = 1.0;
-    setCoverages(cov.data());
-    m_logsize.resize(m_kk);
-    for (size_t k = 0; k < m_kk; k++) {
-        m_logsize[k] = log(size(k));
+    bool added = ThermoPhase::addSpecies(spec);
+    if (added) {
+        m_h0.push_back(0.0);
+        m_s0.push_back(0.0);
+        m_cp0.push_back(0.0);
+        m_mu0.push_back(0.0);
+        m_work.push_back(0.0);
+        m_logsize.push_back(log(size(m_kk-1)));
+        if (m_kk == 1) {
+            vector_fp cov{1.0};
+            setCoverages(cov.data());
+        }
     }
+    return added;
 }
 
 void SurfPhase::setSiteDensity(doublereal n0)
diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp
index 6ce1f801a..a812c1125 100644
--- a/src/thermo/ThermoPhase.cpp
+++ b/src/thermo/ThermoPhase.cpp
@@ -689,17 +689,11 @@ void ThermoPhase::getReferenceComposition(doublereal* const x) const
 
 void ThermoPhase::initThermo()
 {
-    // Check to see that there is at least one species defined in the phase
-    if (m_kk == 0) {
-        throw CanteraError("ThermoPhase::initThermo()",
-                           "Number of species is equal to zero");
-    }
     // Check to see that all of the species thermo objects have been initialized
     if (!m_spthermo->ready(m_kk)) {
         throw CanteraError("ThermoPhase::initThermo()",
                            "Missing species thermo data");
     }
-    xMol_Ref.resize(m_kk, 0.0);
 }
 void ThermoPhase::installSlavePhases(XML_Node* phaseNode)
 {
@@ -711,6 +705,7 @@ bool ThermoPhase::addSpecies(shared_ptr<Species> spec)
     if (added) {
         spec->thermo->validate(spec->name);
         m_spthermo->install_STIT(m_kk-1, spec->thermo);
+        xMol_Ref.push_back(0.0);
     }
     return added;
 }