diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx
index d28ee47ce0..f225f078d3 100644
--- a/rmgpy/chemkin.pyx
+++ b/rmgpy/chemkin.pyx
@@ -1867,7 +1867,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals
         for species, cov_params in kinetics.coverage_dependence.items():
             label = get_species_identifier(species)
             string += f'    COV / {label:<41} '
-            string += f"{cov_params['a']:<9.3g} {cov_params['m']:<9.3g} {cov_params['E'].value_si/4184.:<9.3f} /\n"
+            string += f"{cov_params['a'].value_si:<9.3g} {cov_params['m'].value_si:<9.3g} {cov_params['E'].value_si/4184.:<9.3f} /\n"
 
     if isinstance(kinetics, (_kinetics.ThirdBody, _kinetics.Lindemann, _kinetics.Troe)):
         # Write collider efficiencies
diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index e04cd82e2e..9d0341d2a6 100644
--- a/rmgpy/rmg/input.py
+++ b/rmgpy/rmg/input.py
@@ -105,7 +105,8 @@ def database(
 def catalyst_properties(bindingEnergies=None,
                         surfaceSiteDensity=None,
                         metal=None,
-                        coverageDependence=False):
+                        coverageDependence=False,
+                        thermoCoverageDependence=False,):
     """
     Specify the properties of the catalyst.
     Binding energies of C,H,O,N atoms, and the surface site density.
@@ -152,6 +153,7 @@ def catalyst_properties(bindingEnergies=None,
     else:
         logging.info("Coverage dependence is turned OFF")
     rmg.coverage_dependence = coverageDependence
+    rmg.thermo_coverage_dependence = thermoCoverageDependence
 
 def convert_binding_energies(binding_energies):
     """
@@ -1062,7 +1064,8 @@ def surface_reactor(temperature,
                             sensitive_species=sensitive_species,
                             sensitivity_threshold=sensitivityThreshold,
                             sens_conditions=sens_conditions,
-                            coverage_dependence=rmg.coverage_dependence)
+                            coverage_dependence=rmg.coverage_dependence,
+                            thermo_coverage_dependence=rmg.thermo_coverage_dependence)
     rmg.reaction_systems.append(system)
     system.log_initial_conditions(number=len(rmg.reaction_systems))
 
diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py
index 8df771980d..1e575bf800 100644
--- a/rmgpy/rmg/main.py
+++ b/rmgpy/rmg/main.py
@@ -50,6 +50,7 @@
 import yaml
 from cantera import ck2yaml
 from scipy.optimize import brute
+import cantera as ct
 
 import rmgpy.util as util
 from rmgpy.rmg.model import Species, CoreEdgeReactionModel
@@ -190,6 +191,7 @@ def clear(self):
         self.surface_site_density = None
         self.binding_energies = None
         self.coverage_dependence = False
+        self.thermo_coverage_dependence = False
         self.forbidden_structures = []
 
         self.reaction_model = None
@@ -274,6 +276,7 @@ def load_input(self, path=None):
             self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
             self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si
         self.reaction_model.coverage_dependence = self.coverage_dependence
+        self.reaction_model.thermo_coverage_dependence = self.thermo_coverage_dependence
 
         self.reaction_model.verbose_comments = self.verbose_comments
         self.reaction_model.save_edge_species = self.save_edge_species
@@ -1179,6 +1182,44 @@ def execute(self, initialize=True, **kwargs):
                     os.path.join(self.output_directory, "chemkin", "chem_annotated-gas.inp"),
                     surface_file=(os.path.join(self.output_directory, "chemkin", "chem_annotated-surface.inp")),
                 )
+                if self.thermo_coverage_dependence:
+                    # add thermo coverage dependence to Cantera files
+                    chem_yaml_path = os.path.join(self.output_directory, "cantera", "chem.yaml")
+                    gas = ct.Solution(chem_yaml_path, "gas")
+                    surf = ct.Interface(chem_yaml_path, "surface1", [gas])
+                    with open(chem_yaml_path, 'r') as f:
+                        content = yaml.load(f, Loader=yaml.FullLoader)
+                    
+                    content['phases'][1]['reference-state-coverage'] = 0.11
+                    content['phases'][1]['thermo'] = 'coverage-dependent-surface'
+                    
+                    for s in self.reaction_model.core.species:
+                        if s.contains_surface_site() and s.thermo.thermo_coverage_dependence:
+                            for dep_sp, parameters in s.thermo.thermo_coverage_dependence.items():
+                                mol = Molecule().from_adjacency_list(dep_sp)
+                                for sp in self.reaction_model.core.species:
+                                    if sp.is_isomorphic(mol, strict=False):
+                                        parameters['units'] = {'energy':'J', 'quantity':'mol'}
+                                        parameters['enthalpy-coefficients'] = [value.value_si for value in parameters['enthalpy-coefficients']]
+                                        parameters['entropy-coefficients'] = [value.value_si for value in parameters['entropy-coefficients']]
+                                        try:
+                                            content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'][sp.to_chemkin()] = parameters
+                                        except KeyError:
+                                            content["species"][gas.n_species+surf.species_index(sp.to_chemkin())]['coverage-dependencies'] = {sp.to_chemkin(): parameters}
+
+                    annotated_yaml_path = os.path.join(self.output_directory, "cantera", "chem_annotated.yaml")
+                    with open(annotated_yaml_path, 'r') as f:
+                        annotated_content = yaml.load(f, Loader=yaml.FullLoader)
+                    
+                    annotated_content['phases'] = content['phases']
+                    annotated_content['species'] = content['species']
+
+                    with open(chem_yaml_path, 'w') as output_f:
+                        yaml.dump(content, output_f, sort_keys=False)
+                    
+                    with open(annotated_yaml_path, 'w') as output_f:
+                        yaml.dump(annotated_content, output_f, sort_keys=False)
+                                       
             else:  # gas phase only
                 self.generate_cantera_files(os.path.join(self.output_directory, "chemkin", "chem.inp"))
                 self.generate_cantera_files(os.path.join(self.output_directory, "chemkin", "chem_annotated.inp"))
diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx
index 105a1d68f7..2fad5164f1 100644
--- a/rmgpy/solver/surface.pyx
+++ b/rmgpy/solver/surface.pyx
@@ -43,6 +43,8 @@ cimport rmgpy.constants as constants
 from rmgpy.quantity import Quantity
 from rmgpy.quantity cimport ScalarQuantity
 from rmgpy.solver.base cimport ReactionSystem
+import copy
+from rmgpy.molecule import Molecule
 
 cdef class SurfaceReactor(ReactionSystem):
     """
@@ -66,9 +68,13 @@ cdef class SurfaceReactor(ReactionSystem):
     cdef public ScalarQuantity surface_site_density
     cdef public np.ndarray reactions_on_surface  # (catalyst surface, not core/edge surface)
     cdef public np.ndarray species_on_surface  # (catalyst surface, not core/edge surface)
+    cdef public np.ndarray thermo_coeff_matrix
+    cdef public np.ndarray stoi_matrix
 
     cdef public bint coverage_dependence
     cdef public dict coverage_dependencies
+    cdef public bint thermo_coverage_dependence
+
 
 
     def __init__(self,
@@ -84,6 +90,7 @@ cdef class SurfaceReactor(ReactionSystem):
                  sensitivity_threshold=1e-3,
                  sens_conditions=None,
                  coverage_dependence=False,
+                 thermo_coverage_dependence=False,
                  ):
         ReactionSystem.__init__(self,
                                 termination,
@@ -103,6 +110,7 @@ cdef class SurfaceReactor(ReactionSystem):
         self.surface_volume_ratio = Quantity(surface_volume_ratio)
         self.surface_site_density = Quantity(surface_site_density)
         self.coverage_dependence = coverage_dependence
+        self.thermo_coverage_dependence = thermo_coverage_dependence
         self.V = 0  # will be set from ideal gas law in initialize_model
         self.constant_volume = True
         self.sens_conditions = sens_conditions
@@ -166,6 +174,10 @@ cdef class SurfaceReactor(ReactionSystem):
                                        )
         cdef np.ndarray[np.int_t, ndim=1] species_on_surface, reactions_on_surface
         cdef Py_ssize_t index
+        cdef np.ndarray thermo_coeff_matrix = np.zeros((len(self.species_index), len(self.species_index), 6), dtype=np.float64)
+        cdef np.ndarray stoi_matrix = np.zeros((self.reactant_indices.shape[0], len(self.species_index)), dtype=np.float64)
+        if self.thermo_coverage_dependence:
+            self.thermo_coeff_matrix = thermo_coeff_matrix
         #: 1 if it's on a surface, 0 if it's in the gas phase
         reactions_on_surface = np.zeros((self.num_core_reactions + self.num_edge_reactions), int)
         species_on_surface = np.zeros((self.num_core_species), int)
@@ -195,6 +207,49 @@ cdef class SurfaceReactor(ReactionSystem):
                 means that Species with index 2 in the current simulation is used in
                 Reaction 3 with parameters a=0.1, m=-1, E=12 kJ/mol
                 """
+        for sp, sp_index in self.species_index.items():
+            if sp.contains_surface_site():
+                if self.thermo_coverage_dependence and sp.thermo.thermo_coverage_dependence:
+                    for spec, parameters in sp.thermo.thermo_coverage_dependence.items():
+                        molecule = Molecule().from_adjacency_list(spec)
+                        for species in self.species_index.keys():
+                            if species.is_isomorphic(molecule, strict=False):
+                                species_index = self.species_index[species]
+                                enthalpy_coeff = np.array([p.value_si for p in parameters['enthalpy-coefficients']])
+                                entropy_coeff = np.array([p.value_si for p in parameters['entropy-coefficients']])
+                                thermo_polynomials = np.concatenate((enthalpy_coeff, entropy_coeff), axis=0)
+                                self.thermo_coeff_matrix[sp_index, species_index] = [x for x in thermo_polynomials]
+        # create a stoichiometry matrix for reaction enthalpy and entropy correction 
+        # due to thermodynamic coverage dependence
+        if self.thermo_coverage_dependence:
+            ir = self.reactant_indices
+            ip = self.product_indices
+            for rxn_id, rxn_stoi_num in enumerate(stoi_matrix):
+                if ir[rxn_id, 0] >= self.num_core_species or ir[rxn_id, 1] >= self.num_core_species or ir[rxn_id, 2] >= self.num_core_species:
+                    continue
+                elif ip[rxn_id, 0] >= self.num_core_species or ip[rxn_id, 1] >= self.num_core_species or ip[rxn_id, 2] >= self.num_core_species:
+                    continue
+                else:
+                    if ir[rxn_id, 1] == -1:  # only one reactant
+                        rxn_stoi_num[ir[rxn_id, 0]] += -1
+                    elif ir[rxn_id, 2] == -1:  # only two reactants
+                        rxn_stoi_num[ir[rxn_id, 0]] += -1
+                        rxn_stoi_num[ir[rxn_id, 1]] += -1
+                    else:  # three reactants
+                        rxn_stoi_num[ir[rxn_id, 0]] += -1
+                        rxn_stoi_num[ir[rxn_id, 1]] += -1
+                        rxn_stoi_num[ir[rxn_id, 2]] += -1
+                    if ip[rxn_id, 1] == -1:  # only one product
+                        rxn_stoi_num[ip[rxn_id, 0]] += 1
+                    elif ip[rxn_id, 2] == -1:  # only two products
+                        rxn_stoi_num[ip[rxn_id, 0]] += 1
+                        rxn_stoi_num[ip[rxn_id, 1]] += 1
+                    else:  # three products
+                        rxn_stoi_num[ip[rxn_id, 0]] += 1
+                        rxn_stoi_num[ip[rxn_id, 1]] += 1
+                        rxn_stoi_num[ip[rxn_id, 2]] += 1
+            self.stoi_matrix = stoi_matrix
+        
         self.species_on_surface = species_on_surface
         self.reactions_on_surface = reactions_on_surface
 
@@ -378,9 +433,6 @@ cdef class SurfaceReactor(ReactionSystem):
         cdef np.ndarray[np.float64_t, ndim=2] jacobian, dgdk
         cdef list list_of_coverage_deps
         cdef double surface_site_fraction, total_sites, a, m, E
-
-
-
         ir = self.reactant_indices
         ip = self.product_indices
         equilibrium_constants = self.Keq
@@ -414,7 +466,6 @@ cdef class SurfaceReactor(ReactionSystem):
         V = self.V  # constant volume reactor
         A = self.V * surface_volume_ratio_si  # area
         total_sites = self.surface_site_density.value_si * A  # todo: double check units
-
         for j in range(num_core_species):
             if species_on_surface[j]:
                 C[j] = (N[j] / V) / surface_volume_ratio_si
@@ -422,6 +473,26 @@ cdef class SurfaceReactor(ReactionSystem):
                 C[j] = N[j] / V
             #: surface species are in mol/m2, gas phase are in mol/m3
             core_species_concentrations[j] = C[j]
+        
+        # Thermodynamic coverage dependence
+        if self.thermo_coverage_dependence:
+            coverages = []
+            for i in range(len(N)):
+                if species_on_surface[i]:
+                    surface_site_fraction = N[i] / total_sites
+                else:
+                    surface_site_fraction = 0
+                coverages.append(surface_site_fraction)
+            coverages = np.array(coverages)
+            thermo_dep_coverage = np.stack([coverages, coverages**2, coverages**3, -self.T.value_si*coverages, -self.T.value_si*coverages**2, -self.T.value_si*coverages**3])
+            free_energy_coverage_corrections = []
+            for matrix in self.thermo_coeff_matrix:
+                sp_free_energy_correction = np.diag(np.dot(matrix, thermo_dep_coverage)).sum()
+                free_energy_coverage_corrections.append(sp_free_energy_correction)
+            rxns_free_energy_change = np.diag(np.dot(self.stoi_matrix, np.transpose(np.array([free_energy_coverage_corrections]))))
+            corrected_K_eq = copy.deepcopy(self.Keq)
+            corrected_K_eq *= np.exp(-1 * rxns_free_energy_change / (constants.R * self.T.value_si))
+            kr = kf / corrected_K_eq
 
         # Coverage dependence
         coverage_corrections = np.ones_like(kf, float)
diff --git a/rmgpy/thermo/nasa.pxd b/rmgpy/thermo/nasa.pxd
index b462301f72..6922c51e54 100644
--- a/rmgpy/thermo/nasa.pxd
+++ b/rmgpy/thermo/nasa.pxd
@@ -56,6 +56,7 @@ cdef class NASAPolynomial(HeatCapacityModel):
 cdef class NASA(HeatCapacityModel):
 
     cdef public NASAPolynomial poly1, poly2, poly3
+    cdef public dict _thermo_coverage_dependence
     
     cpdef NASAPolynomial select_polynomial(self, double T)
 
diff --git a/rmgpy/thermo/nasa.pyx b/rmgpy/thermo/nasa.pyx
index e484f3036f..0b04438654 100644
--- a/rmgpy/thermo/nasa.pyx
+++ b/rmgpy/thermo/nasa.pyx
@@ -34,6 +34,8 @@ cimport numpy as np
 from libc.math cimport log
 
 cimport rmgpy.constants as constants
+from rmgpy.util import np_list
+import rmgpy.quantity as quantity
 
 ################################################################################
 
@@ -43,15 +45,15 @@ cdef class NASAPolynomial(HeatCapacityModel):
     seven-coefficient and nine-coefficient variations are supported.
     The attributes are:
     
-    =============== ============================================================
-    Attribute       Description
-    =============== ============================================================
-    `coeffs`        The seven or nine NASA polynomial coefficients
-    `Tmin`          The minimum temperature in K at which the model is valid, or zero if unknown or undefined
-    `Tmax`          The maximum temperature in K at which the model is valid, or zero if unknown or undefined
-    `E0`            The energy at zero Kelvin (including zero point energy)
-    `comment`       Information about the model (e.g. its source)
-    =============== ============================================================
+    =========  ============================================================
+    Attribute  Description
+    =========  ============================================================
+    `coeffs`   The seven or nine NASA polynomial coefficients
+    `Tmin`     The minimum temperature in K at which the model is valid, or zero if unknown or undefined
+    `Tmax`     The maximum temperature in K at which the model is valid, or zero if unknown or undefined
+    `E0`       The energy at zero Kelvin (including zero point energy)
+    `comment`  Information about the model (e.g. its source)
+    =========  ============================================================
 
     """
     
@@ -203,24 +205,26 @@ cdef class NASA(HeatCapacityModel):
     A heat capacity model based on a set of one, two, or three 
     :class:`NASAPolynomial` objects. The attributes are:
     
-    =============== ============================================================
-    Attribute       Description
-    =============== ============================================================
-    `polynomials`   The list of NASA polynomials to use in this model
-    `Tmin`          The minimum temperature in K at which the model is valid, or zero if unknown or undefined
-    `Tmax`          The maximum temperature in K at which the model is valid, or zero if unknown or undefined
-    `E0`            The energy at zero Kelvin (including zero point energy)
-    `comment`       Information about the model (e.g. its source)
-    =============== ============================================================
+    ============================ ============================================================
+    Attribute                    Description
+    ============================ ============================================================
+    `polynomials`                The list of NASA polynomials to use in this model
+    `Tmin`                       The minimum temperature in K at which the model is valid, or zero if unknown or undefined
+    `Tmax`                       The maximum temperature in K at which the model is valid, or zero if unknown or undefined
+    `E0`                         The energy at zero Kelvin (including zero point energy)
+    `thermo_coverage_dependence` The coverage dependence of the thermo
+    `comment`                    Information about the model (e.g. its source)
+    ============================ ============================================================
 
     """
     
-    def __init__(self, polynomials=None, Tmin=None, Tmax=None, E0=None, Cp0=None, CpInf=None, label='', comment=''):
+    def __init__(self, polynomials=None, Tmin=None, Tmax=None, E0=None, Cp0=None, CpInf=None, label='', thermo_coverage_dependence=None, comment=''):
         HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, E0=E0, Cp0=Cp0, CpInf=CpInf, label=label, comment=comment)
         if isinstance(polynomials, dict):
             self.polynomials = list(polynomials.values())
         else:
             self.polynomials = polynomials
+        self.thermo_coverage_dependence = thermo_coverage_dependence
     
     def __repr__(self):
         """
@@ -235,6 +239,7 @@ cdef class NASA(HeatCapacityModel):
         if self.Cp0 is not None: string += ', Cp0={0!r}'.format(self.Cp0)
         if self.CpInf is not None: string += ', CpInf={0!r}'.format(self.CpInf)
         if self.label != '': string += ', label="""{0}"""'.format(self.label)
+        if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence)
         if self.comment != '': string += ', comment="""{0}"""'.format(self.comment)
         string += ')'
         return string
@@ -243,7 +248,7 @@ cdef class NASA(HeatCapacityModel):
         """
         A helper function used when pickling an object.
         """
-        return (NASA, (self.polynomials, self.Tmin, self.Tmax, self.E0, self.Cp0, self.CpInf, self.label, self.comment))
+        return (NASA, (self.polynomials, self.Tmin, self.Tmax, self.E0, self.Cp0, self.CpInf, self.label, self.thermo_coverage_dependence, self.comment))
 
     cpdef dict as_dict(self):
         """
@@ -300,6 +305,21 @@ cdef class NASA(HeatCapacityModel):
         else:
             raise ValueError('No valid NASA polynomial at temperature {0:g} K.'.format(T))
     
+    property thermo_coverage_dependence:
+        """The coverage dependence of the thermo"""
+        def __get__(self):
+            return self._thermo_coverage_dependence
+        def __set__(self, value):
+            self._thermo_coverage_dependence = {}
+            if value:
+                 for species, parameters in value.items():
+                    # just the polynomial model for now
+                     processed_parameters = {'model': parameters['model'],
+                                             'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]),
+                                             'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]),
+                                             }
+                     self._thermo_coverage_dependence[species] = processed_parameters
+    
     cpdef double get_heat_capacity(self, double T) except -1000000000:
         """
         Return the constant-pressure heat capacity
@@ -347,6 +367,7 @@ cdef class NASA(HeatCapacityModel):
             Cp0 = self.Cp0,
             CpInf = self.CpInf,
             E0 = self.E0,
+            thermo_coverage_dependence = self.thermo_coverage_dependence,
             comment = self.comment
         )
 
@@ -379,7 +400,7 @@ cdef class NASA(HeatCapacityModel):
         H298 = self.get_enthalpy(298)
         S298 = self.get_entropy(298)
         
-        return Wilhoit(label=self.label,comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298)
+        return Wilhoit(label=self.label, thermo_coverage_dependence=self.thermo_coverage_dependence, comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298)
     
     cpdef NASA change_base_enthalpy(self, double deltaH):
         """
@@ -399,6 +420,7 @@ cdef class NASA(HeatCapacityModel):
             poly.change_base_entropy(deltaS)
         return self
 
+    # need to modify this to include the thermo coverage dependence
     def to_cantera(self):
         """
         Return the cantera equivalent NasaPoly2 object from this NASA object.
diff --git a/rmgpy/thermo/thermodata.pxd b/rmgpy/thermo/thermodata.pxd
index d9052db0b4..0a1d08b12f 100644
--- a/rmgpy/thermo/thermodata.pxd
+++ b/rmgpy/thermo/thermodata.pxd
@@ -38,6 +38,7 @@ cdef class ThermoData(HeatCapacityModel):
     
     cdef public ScalarQuantity _H298, _S298
     cdef public ArrayQuantity _Tdata, _Cpdata
+    cdef public dict _thermo_coverage_dependence
     
     cpdef double get_heat_capacity(self, double T) except -1000000000
 
diff --git a/rmgpy/thermo/thermodata.pyx b/rmgpy/thermo/thermodata.pyx
index 5248e79a56..6797fa99e4 100644
--- a/rmgpy/thermo/thermodata.pyx
+++ b/rmgpy/thermo/thermodata.pyx
@@ -35,6 +35,7 @@ cimport numpy as np
 from libc.math cimport log
 
 import rmgpy.quantity as quantity
+from rmgpy.util import np_list
 
 ################################################################################
 
@@ -43,27 +44,29 @@ cdef class ThermoData(HeatCapacityModel):
     A heat capacity model based on a set of discrete heat capacity data points.
     The attributes are:
     
-    =============== ============================================================
-    Attribute       Description
-    =============== ============================================================
-    `Tdata`         An array of temperatures at which the heat capacity is known
-    `Cpdata`        An array of heat capacities at the given temperatures
-    `H298`          The standard enthalpy of formation at 298 K
-    `S298`          The standard entropy at 298 K
-    `Tmin`          The minimum temperature at which the model is valid, or zero if unknown or undefined
-    `Tmax`          The maximum temperature at which the model is valid, or zero if unknown or undefined
-    `E0`            The energy at zero Kelvin (including zero point energy)
-    `comment`       Information about the model (e.g. its source)
-    =============== ============================================================
+    ===========================  ============================================================
+    Attribute                    Description
+    ===========================  ============================================================
+    `Tdata`                      An array of temperatures at which the heat capacity is known
+    `Cpdata`                     An array of heat capacities at the given temperatures
+    `H298`                       The standard enthalpy of formation at 298 K
+    `S298`                       The standard entropy at 298 K
+    `Tmin`                       The minimum temperature at which the model is valid, or zero if unknown or undefined
+    `Tmax`                       The maximum temperature at which the model is valid, or zero if unknown or undefined
+    `E0`                         The energy at zero Kelvin (including zero point energy)
+    `thermo_coverage_dependence` The coverage dependence of the thermo
+    `comment`                    Information about the model (e.g. its source)
+    ===========================  ============================================================
     
     """
 
-    def __init__(self, Tdata=None, Cpdata=None, H298=None, S298=None, Cp0=None, CpInf=None, Tmin=None, Tmax=None, E0=None, label = '',comment=''):
+    def __init__(self, Tdata=None, Cpdata=None, H298=None, S298=None, Cp0=None, CpInf=None, Tmin=None, Tmax=None, E0=None, label = '', thermo_coverage_dependence=None, comment=''):
         HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, E0=E0, Cp0=Cp0, CpInf=CpInf, label=label, comment=comment)
         self.H298 = H298
         self.S298 = S298
         self.Tdata = Tdata
         self.Cpdata = Cpdata
+        self.thermo_coverage_dependence = thermo_coverage_dependence
     
     def __repr__(self):
         """
@@ -77,6 +80,7 @@ cdef class ThermoData(HeatCapacityModel):
         if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax)
         if self.E0 is not None: string += ', E0={0!r}'.format(self.E0)
         if self.label != '': string +=', label="""{0}"""'.format(self.label)
+        if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence)
         if self.comment != '': string += ', comment="""{0}"""'.format(self.comment)
         string += ')'
         return string
@@ -85,7 +89,7 @@ cdef class ThermoData(HeatCapacityModel):
         """
         A helper function used when pickling a ThermoData object.
         """
-        return (ThermoData, (self.Tdata, self.Cpdata, self.H298, self.S298, self.Cp0, self.CpInf, self.Tmin, self.Tmax, self.E0, self.label, self.comment))
+        return (ThermoData, (self.Tdata, self.Cpdata, self.H298, self.S298, self.Cp0, self.CpInf, self.Tmin, self.Tmax, self.E0, self.label, self.thermo_coverage_dependence, self.comment))
 
     property Tdata:
         """An array of temperatures at which the heat capacity is known."""
@@ -114,6 +118,21 @@ cdef class ThermoData(HeatCapacityModel):
             return self._S298
         def __set__(self, value):
             self._S298 = quantity.Entropy(value)
+    
+    property thermo_coverage_dependence:
+        """The coverage dependence of the thermo"""
+        def __get__(self):
+            return self._thermo_coverage_dependence
+        def __set__(self, value):
+            self._thermo_coverage_dependence = {}
+            if value:
+                 for species, parameters in value.items():
+                    # just the polynomial model for now
+                     processed_parameters = {'model': parameters['model'],
+                                             'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]),
+                                             'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]),
+                                             }
+                     self._thermo_coverage_dependence[species] = processed_parameters
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -361,12 +380,12 @@ cdef class ThermoData(HeatCapacityModel):
         S298 = self.get_entropy(298)
         Cp0 = self._Cp0.value_si
         CpInf = self._CpInf.value_si
-        
+        thermo_cov_dep = self.thermo_coverage_dependence
         
         if B:
-            return Wilhoit(label=self.label,comment=self.comment).fit_to_data_for_constant_b(Tdata, Cpdata, Cp0, CpInf, H298, S298, B=B)
+            return Wilhoit(label=self.label, thermo_coverage_dependence=thermo_cov_dep, comment=self.comment).fit_to_data_for_constant_b(Tdata, Cpdata, Cp0, CpInf, H298, S298, B=B)
         else:
-            return Wilhoit(label=self.label,comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298)
+            return Wilhoit(label=self.label, thermo_coverage_dependence=thermo_cov_dep, comment=self.comment).fit_to_data(Tdata, Cpdata, Cp0, CpInf, H298, S298)
 
     cpdef NASA to_nasa(self, double Tmin, double Tmax, double Tint, bint fixedTint=False, bint weighting=True, int continuity=3):
         """
diff --git a/rmgpy/thermo/wilhoit.pxd b/rmgpy/thermo/wilhoit.pxd
index 576c8a2e67..e2cd2419e4 100644
--- a/rmgpy/thermo/wilhoit.pxd
+++ b/rmgpy/thermo/wilhoit.pxd
@@ -38,6 +38,7 @@ cdef class Wilhoit(HeatCapacityModel):
     
     cdef public ScalarQuantity _B, _H0, _S0
     cdef public double a0, a1, a2, a3
+    cdef public dict _thermo_coverage_dependence
 
     cpdef dict as_dict(self)
 
diff --git a/rmgpy/thermo/wilhoit.pyx b/rmgpy/thermo/wilhoit.pyx
index 20a1061b18..3342f64c48 100644
--- a/rmgpy/thermo/wilhoit.pyx
+++ b/rmgpy/thermo/wilhoit.pyx
@@ -35,6 +35,7 @@ from libc.math cimport sqrt, log
 
 cimport rmgpy.constants as constants
 import rmgpy.quantity as quantity
+from rmgpy.util import np_list
 
 # Prior to numpy 1.14, `numpy.linalg.lstsq` does not accept None as a value
 RCOND = -1 if int(np.__version__.split('.')[1]) < 14 else None
@@ -45,24 +46,25 @@ cdef class Wilhoit(HeatCapacityModel):
     """
     A heat capacity model based on the Wilhoit equation. The attributes are:
     
-    =============== ============================================================
-    Attribute       Description
-    =============== ============================================================
-    `a0`            The zeroth-order Wilhoit polynomial coefficient
-    `a1`            The first-order Wilhoit polynomial coefficient
-    `a2`            The second-order Wilhoit polynomial coefficient
-    `a3`            The third-order Wilhoit polynomial coefficient
-    `H0`            The integration constant for enthalpy (not H at T=0)
-    `S0`            The integration constant for entropy (not S at T=0)
-    `E0`            The energy at zero Kelvin (including zero point energy)
-    `B`             The Wilhoit scaled temperature coefficient in K
-    `Tmin`          The minimum temperature in K at which the model is valid, or zero if unknown or undefined
-    `Tmax`          The maximum temperature in K at which the model is valid, or zero if unknown or undefined
-    `comment`       Information about the model (e.g. its source)
-    =============== ============================================================
+    ============================= =========================================================================================
+    Attribute                     Description
+    ============================= =========================================================================================
+    `a0`                          The zeroth-order Wilhoit polynomial coefficient
+    `a1`                          The first-order Wilhoit polynomial coefficient
+    `a2`                          The second-order Wilhoit polynomial coefficient
+    `a3`                          The third-order Wilhoit polynomial coefficient
+    `H0`                          The integration constant for enthalpy (not H at T=0)
+    `S0`                          The integration constant for entropy (not S at T=0)
+    `E0`                          The energy at zero Kelvin (including zero point energy)
+    `B`                           The Wilhoit scaled temperature coefficient in K
+    `Tmin`                        The minimum temperature in K at which the model is valid, or zero if unknown or undefined
+    `Tmax`                        The maximum temperature in K at which the model is valid, or zero if unknown or undefined
+    `thermo_covreage_dependence`  The coverage dependence of the thermo
+    `comment`                     Information about the model (e.g. its source)
+    ============================= =========================================================================================
     """
 
-    def __init__(self, Cp0=None, CpInf=None, a0=0.0, a1=0.0, a2=0.0, a3=0.0, H0=None, S0=None, B=None, Tmin=None, Tmax=None, label='', comment=''):
+    def __init__(self, Cp0=None, CpInf=None, a0=0.0, a1=0.0, a2=0.0, a3=0.0, H0=None, S0=None, B=None, Tmin=None, Tmax=None, label='', thermo_coverage_dependence=None, comment=''):
         HeatCapacityModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Cp0=Cp0, CpInf=CpInf, label=label, comment=comment)
         self.B = B
         self.a0 = a0
@@ -71,6 +73,7 @@ cdef class Wilhoit(HeatCapacityModel):
         self.a3 = a3
         self.H0 = H0
         self.S0 = S0
+        self.thermo_coverage_dependence = thermo_coverage_dependence
     
     def __repr__(self):
         """
@@ -82,6 +85,7 @@ cdef class Wilhoit(HeatCapacityModel):
         if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin)
         if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax)
         if self.label != '': string += ', label="""{0}"""'.format(self.label)
+        if self.thermo_coverage_dependence is not None: string += ', thermo_coverage_dependence={0!r}'.format(self.thermo_coverage_dependence)
         if self.comment != '': string += ', comment="""{0}"""'.format(self.comment)
         string += ')'
         return string
@@ -90,7 +94,7 @@ cdef class Wilhoit(HeatCapacityModel):
         """
         A helper function used when pickling a Wilhoit object.
         """
-        return (Wilhoit, (self.Cp0, self.CpInf, self.a0, self.a1, self.a2, self.a3, self.H0, self.S0, self.B, self.Tmin, self.Tmax, self.label, self.comment))
+        return (Wilhoit, (self.Cp0, self.CpInf, self.a0, self.a1, self.a2, self.a3, self.H0, self.S0, self.B, self.Tmin, self.Tmax, self.label, self.thermo_coverage_dependence, self.comment))
 
     cpdef dict as_dict(self):
         """
@@ -136,6 +140,21 @@ cdef class Wilhoit(HeatCapacityModel):
             return self._S0
         def __set__(self, value):
             self._S0 = quantity.Entropy(value)
+    
+    property thermo_coverage_dependence:
+        """The coverage dependence of the thermo"""
+        def __get__(self):
+            return self._thermo_coverage_dependence
+        def __set__(self, value):
+            self._thermo_coverage_dependence = {}
+            if value:
+                 for species, parameters in value.items():
+                    # just the polynomial model for now
+                     processed_parameters = {'model': parameters['model'],
+                                             'enthalpy-coefficients': np_list([quantity.Enthalpy(p) for p in parameters['enthalpy-coefficients']]),
+                                             'entropy-coefficients': np_list([quantity.Entropy(p) for p in parameters['entropy-coefficients']]),
+                                             }
+                     self._thermo_coverage_dependence[species] = processed_parameters
 
     cpdef double get_heat_capacity(self, double T) except -1000000000:
         """
@@ -195,7 +214,8 @@ cdef class Wilhoit(HeatCapacityModel):
             self.Cp0, self.CpInf, 
             self.a0, self.a1, self.a2, self.a3, 
             self.H0, self.S0, self.B, 
-            Tmin=self.Tmin, Tmax=self.Tmax, comment=self.comment,
+            Tmin=self.Tmin, Tmax=self.Tmax, thermo_coverage_dependence=self.thermo_coverage_dependence,
+            comment=self.comment,
         )
     
     @cython.boundscheck(False)
@@ -478,6 +498,7 @@ cdef class Wilhoit(HeatCapacityModel):
             Cp0 = self.Cp0,
             CpInf = self.CpInf,
             E0 = self.E0,
+            thermo_coverage_dependence = self.thermo_coverage_dependence,
             comment = self.comment
         )
     
@@ -582,6 +603,7 @@ cdef class Wilhoit(HeatCapacityModel):
             Cp0 = self.Cp0,
             CpInf = self.CpInf,
             label = self.label,
+            thermo_coverage_dependence = self.thermo_coverage_dependence,
             comment = self.comment,
         )
     
diff --git a/rmgpy/util.py b/rmgpy/util.py
index 15769e02e6..365833379e 100644
--- a/rmgpy/util.py
+++ b/rmgpy/util.py
@@ -33,6 +33,7 @@
 import shutil
 import time
 from functools import wraps
+import numpy as np
 
 
 class Subject(object):
@@ -238,3 +239,13 @@ def as_list(item, default=None):
         return default
     else:
         return [item]
+
+class np_list(np.ndarray):
+    """
+    A subclass of numpy.ndarray which rendered as a list when printed.
+    """
+    def __new__(cls, input_array):
+        obj = np.asarray(input_array).view(cls)
+        return obj
+    def __repr__(self):
+        return str(self.tolist())
diff --git a/test/rmgpy/solver/solverSurfaceTest.py b/test/rmgpy/solver/solverSurfaceTest.py
index 36c08b4ad8..42fb9c8140 100644
--- a/test/rmgpy/solver/solverSurfaceTest.py
+++ b/test/rmgpy/solver/solverSurfaceTest.py
@@ -762,3 +762,416 @@ def test_solve_ch3_coverage_dependence(self):
         # Check that coverages are different
         assert not np.allclose(y, y_off)
         assert not np.allclose(species_rates, species_rates_off)
+
+    def test_solve_h2_thermo_coverage_dependence(self):
+        """
+        Test the surface batch reactor can properly apply thermodynamic coverage dependent parameters
+        with the dissociative adsorption of H2.
+
+        Here we choose a kinetic model consisting of the dissociative adsorption reaction
+        H2 + 2X <=> 2 HX
+        We use a SurfaceArrhenius for the rate expression.
+        """
+        h2 = Species(
+            molecule=[Molecule().from_smiles("[H][H]")],
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=(
+                    [6.955, 6.955, 6.956, 6.961, 7.003, 7.103, 7.502],
+                    "cal/(mol*K)",
+                ),
+                H298=(0, "kcal/mol"),
+                S298=(31.129, "cal/(mol*K)"),
+            ),
+        )
+        
+        x = Species(
+            molecule=[Molecule().from_adjacency_list("1 X u0 p0")],
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "cal/(mol*K)"),
+                H298=(0.0, "kcal/mol"),
+                S298=(0.0, "cal/(mol*K)"),
+            ),
+        )
+        
+        hx = Species(
+            molecule=[Molecule().from_adjacency_list("1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}")],
+            thermo=ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=([1.50, 2.58, 3.40, 4.00, 4.73, 5.13, 5.57], "cal/(mol*K)"),
+                H298=(-11.26, "kcal/mol"),
+                S298=(0.44, "cal/(mol*K)"),
+                thermo_coverage_dependence={"1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}":{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,5,3]},}
+            ),
+        )
+        
+        rxn1 = Reaction(
+            reactants=[h2, x, x],
+            products=[hx, hx],
+            kinetics=SurfaceArrhenius(
+                A=(9.05e18, "cm^5/(mol^2*s)"),
+                n=0.5,
+                Ea=(5.0, "kJ/mol"),
+                T0=(1.0, "K"),
+            ),
+        )
+
+        rxn2 = Reaction(
+            reactants=[h2, x, x],
+            products=[hx, hx],
+            kinetics=SurfaceArrhenius(
+                A=(9.05e-18, "cm^5/(mol^2*s)"),  # 1e36 times slower
+                n=0.5,
+                Ea=(5.0, "kJ/mol"),
+                T0=(1.0, "K"),
+            ),
+        )
+
+        core_species = [h2, x, hx]
+        edge_species = []
+        core_reactions = [rxn1]
+        edge_reactions = []
+
+        T = 600
+        P_initial = 1.0e5
+        rxn_system = SurfaceReactor(
+            T,
+            P_initial,
+            n_sims=1,
+            initial_gas_mole_fractions={h2: 1.0},
+            initial_surface_coverages={x: 1.0},
+            surface_volume_ratio=(1e1, "m^-1"),
+            surface_site_density=(2.72e-9, "mol/cm^2"),
+            coverage_dependence=True,
+            thermo_coverage_dependence=True,
+            termination=[],
+        )
+
+        rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions)
+
+        tlist = np.logspace(-13, -5, 81, dtype=float)
+
+        assert isinstance(hx.thermo.thermo_coverage_dependence, dict)  # check to make sure coverage_dependence is still the correct type
+        for species, parameters in hx.thermo.thermo_coverage_dependence.items():
+            assert isinstance(species, str)  # species should be an ajacency list
+            assert isinstance(parameters, dict)
+            assert parameters["model"] is not None
+            assert parameters["enthalpy-coefficients"] is not None
+            assert parameters["entropy-coefficients"] is not None
+        assert np.array_equal(rxn_system.stoi_matrix, np.array([[-1., -2.,  2.]]))
+        thermo_coeffs = np.array([np.zeros((3,6))]*3)
+        thermo_coeffs[-1][-1] = [1., 2., 3., 1., 5., 3.]
+        assert np.array_equal(rxn_system.thermo_coeff_matrix, thermo_coeffs)
+        # Integrate to get the solution at each time point
+        t = []
+        y = []
+        reaction_rates = []
+        species_rates = []
+        start_time = time.time()
+        for t1 in tlist:
+            rxn_system.advance(t1)
+            t.append(rxn_system.t)
+            # You must make a copy of y because it is overwritten by DASSL at
+            # each call to advance()
+            y.append(rxn_system.y.copy())
+            reaction_rates.append(rxn_system.core_reaction_rates.copy())
+            species_rates.append(rxn_system.core_species_rates.copy())
+        run_time = time.time() - start_time
+        print(f"Simulation took {run_time:.3e} seconds")
+
+        # Convert the solution vectors to np arrays
+        t = np.array(t, float)
+        y = np.array(y, float)
+        reaction_rates = np.array(reaction_rates, float)
+        species_rates = np.array(species_rates, float)
+        total_sites = y[0, 1]
+
+        # Check that we're computing the species fluxes correctly
+        for i in range(t.shape[0]):
+            assert abs(reaction_rates[i, 0] - -1.0 * species_rates[i, 0]) < abs(
+                1e-6 * reaction_rates[i, 0]
+            )
+            assert abs(reaction_rates[i, 0] - -0.5 * species_rates[i, 1]) < abs(
+                1e-6 * reaction_rates[i, 0]
+            )
+            assert abs(reaction_rates[i, 0] - 0.5 * species_rates[i, 2]) < abs(
+                1e-6 * reaction_rates[i, 0]
+            )
+
+        # Check that we've reached equilibrium
+        assert abs(reaction_rates[-1, 0] - 0.0) < 1e-2
+
+    def test_solve_ch3_thermo_coverage_dependence(self):
+        """
+        Test the surface batch reactor can properly apply thermodynamic coverage dependent parameters
+        with the nondissociative adsorption of CH3
+
+        Here we choose a kinetic model consisting of the  adsorption reaction
+        CH3 + X <=>  CH3X
+        We use a sticking coefficient for the rate expression.
+        """
+
+        ch3 = Species(
+            molecule=[Molecule().from_smiles("[CH3]")],
+            thermo=NASA(
+                polynomials=[
+                    NASAPolynomial(
+                        coeffs=[
+                            3.91547,
+                            0.00184155,
+                            3.48741e-06,
+                            -3.32746e-09,
+                            8.49953e-13,
+                            16285.6,
+                            0.351743,
+                        ],
+                        Tmin=(100, "K"),
+                        Tmax=(1337.63, "K"),
+                    ),
+                    NASAPolynomial(
+                        coeffs=[
+                            3.54146,
+                            0.00476786,
+                            -1.82148e-06,
+                            3.28876e-10,
+                            -2.22545e-14,
+                            16224,
+                            1.66032,
+                        ],
+                        Tmin=(1337.63, "K"),
+                        Tmax=(5000, "K"),
+                    ),
+                ],
+                Tmin=(100, "K"),
+                Tmax=(5000, "K"),
+                E0=(135.382, "kJ/mol"),
+                comment="""Thermo library: primaryThermoLibrary + radical(CH3)""",
+            ),
+            molecular_weight=(15.0345, "amu"),
+        )
+
+        x = Species(
+            molecule=[Molecule().from_adjacency_list("1 X u0 p0")],
+            thermo=NASA(
+                polynomials=[
+                    NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(298, "K"), Tmax=(1000, "K")),
+                    NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(1000, "K"), Tmax=(2000, "K")),
+                ],
+                Tmin=(298, "K"),
+                Tmax=(2000, "K"),
+                E0=(-6.19426, "kJ/mol"),
+                comment="""Thermo library: surfaceThermo""",
+            ),
+        )
+
+        ch3x = Species(
+            molecule=[
+                Molecule().from_adjacency_list(
+                    """1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+4 H u0 p0 c0 {1,S}
+5 X u0 p0 c0 {1,S}"""
+                )
+            ],
+            thermo=NASA(
+                polynomials=[
+                    NASAPolynomial(
+                        coeffs=[
+                            -0.552219,
+                            0.026442,
+                            -3.55617e-05,
+                            2.60044e-08,
+                            -7.52707e-12,
+                            -4433.47,
+                            0.692144,
+                        ],
+                        Tmin=(298, "K"),
+                        Tmax=(1000, "K"),
+                    ),
+                    NASAPolynomial(
+                        coeffs=[
+                            3.62557,
+                            0.00739512,
+                            -2.43797e-06,
+                            1.86159e-10,
+                            3.6485e-14,
+                            -5187.22,
+                            -18.9668,
+                        ],
+                        Tmin=(1000, "K"),
+                        Tmax=(2000, "K"),
+                    ),
+                ],
+                Tmin=(298, "K"),
+                Tmax=(2000, "K"),
+                E0=(-39.1285, "kJ/mol"),
+                thermo_coverage_dependence={"1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} \n 2 H u0 p0 c0 {1,S} \n 3 H u0 p0 c0 {1,S} \n 4 H u0 p0 c0 {1,S} \n 5 X u0 p0 c0 {1,S}":
+                                            {'model':'polynomial', 'enthalpy-coefficients':[1e5,2,3], "entropy-coefficients":[1,5,3]},},
+                comment="""Thermo library: surfaceThermoNi111""",
+            ),
+        )
+
+        rxn1 = Reaction(
+            reactants=[ch3, x],
+            products=[ch3x],
+            kinetics=StickingCoefficient(
+                A=0.1,
+                n=0,
+                Ea=(0, "kcal/mol"),
+                T0=(1, "K"),
+                Tmin=(200, "K"),
+                Tmax=(3000, "K"),
+                comment="""Exact match found for rate rule (Adsorbate;VacantSite)""",
+            ),
+        )
+        core_species = [ch3, x, ch3x]
+        edge_species = []
+        core_reactions = [rxn1]
+        edge_reactions = []
+
+        T = 800.0
+        P_initial = 1.0e5
+        rxn_system = SurfaceReactor(
+            T,
+            P_initial,
+            n_sims=1,
+            initial_gas_mole_fractions={ch3: 1.0},
+            initial_surface_coverages={x: 1.0},
+            surface_volume_ratio=(1.0, "m^-1"),
+            surface_site_density=(2.72e-9, "mol/cm^2"),
+            thermo_coverage_dependence=True,
+            termination=[],
+        )
+        # in chemkin, the sites are mostly occupied in about 1e-8 seconds.
+
+        rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions)
+
+        tlist = np.logspace(-13, -5, 81, dtype=float)
+
+        print("Surface site density:", rxn_system.surface_site_density.value_si)
+
+        print(
+            "rxn1 rate coefficient",
+            rxn1.get_surface_rate_coefficient(rxn_system.T.value_si, rxn_system.surface_site_density.value_si),
+        )
+
+        assert isinstance(ch3x.thermo.thermo_coverage_dependence, dict)  # check to make sure coverage_dependence is still the correct type
+        for species, parameters in ch3x.thermo.thermo_coverage_dependence.items():
+            assert isinstance(species, str)  # species should be an ajacency list
+            assert isinstance(parameters, dict)
+            assert parameters["model"] is not None
+            assert parameters["enthalpy-coefficients"] is not None
+            assert parameters["entropy-coefficients"] is not None
+        
+        # check thermo_coverage_dependence is on
+        # and thermo_coeff_matrix and stoi_matrix are correctly created
+        assert rxn_system.thermo_coverage_dependence is True
+        assert np.array_equal(rxn_system.stoi_matrix, np.array([[-1, -1, 1]], dtype=float))
+        thermo_coeff_matrix = np.array([np.zeros((3,6))]*3)
+        thermo_coeff_matrix[-1][-1] = [1e5, 2, 3, 1, 5, 3]
+        assert np.array_equal(rxn_system.thermo_coeff_matrix, thermo_coeff_matrix)
+
+        # Integrate to get the solution at each time point
+        t = []
+        y = []
+        reaction_rates = []
+        species_rates = []
+        t.append(rxn_system.t)
+        # You must make a copy of y because it is overwritten by DASSL at
+        # each call to advance()
+        y.append(rxn_system.y.copy())
+        reaction_rates.append(rxn_system.core_reaction_rates.copy())
+        species_rates.append(rxn_system.core_species_rates.copy())
+        print("time: ", t)
+        print("moles:", y)
+        print("reaction rates:", reaction_rates)
+        print("species rates:", species_rates)
+        start_time = time.time()
+        for t1 in tlist:
+            rxn_system.advance(t1)
+            t.append(rxn_system.t)
+            # You must make a copy of y because it is overwritten by DASSL at
+            # each call to advance()
+            y.append(rxn_system.y.copy())
+            reaction_rates.append(rxn_system.core_reaction_rates.copy())
+            species_rates.append(rxn_system.core_species_rates.copy())
+        run_time = time.time() - start_time
+        print(f"Simulation took {run_time:.3e} seconds")
+
+        # Convert the solution vectors to np arrays
+        t = np.array(t, float)
+        y = np.array(y, float)
+        reaction_rates = np.array(reaction_rates, float)
+        species_rates = np.array(species_rates, float)
+        V = constants.R * rxn_system.T.value_si * np.sum(y) / rxn_system.P_initial.value_si
+
+        # Check that we're computing the species fluxes correctly
+        for i in range(t.shape[0]):
+            assert abs(reaction_rates[i, 0] - -species_rates[i, 0]) < abs(
+                1e-6 * reaction_rates[i, 0]
+            )
+            assert abs(reaction_rates[i, 0] - -species_rates[i, 1]) < abs(
+                1e-6 * reaction_rates[i, 0]
+            )
+            assert abs(reaction_rates[i, 0] - species_rates[i, 2]) < abs(
+                1e-6 * reaction_rates[i, 0]
+            )
+
+        # Check that we've reached equilibrium by the end
+        assert abs(reaction_rates[-1, 0] - 0.0) < 1e-2
+
+        # Run model with Covdep off so we can test that it is actually being implemented
+        rxn_system = SurfaceReactor(
+            T,
+            P_initial,
+            n_sims=1,
+            initial_gas_mole_fractions={ch3: 1.0},
+            initial_surface_coverages={x: 1.0},
+            surface_volume_ratio=(1.0, "m^-1"),
+            surface_site_density=(2.72e-9, "mol/cm^2"),
+            termination=[],
+        )
+
+        rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions)
+        
+        # check thermo_coverage_dependence is off
+        # and thermo_coeff_matrix and stoi_matrix are not created
+        assert rxn_system.thermo_coverage_dependence is False
+        assert rxn_system.thermo_coeff_matrix is None
+        assert rxn_system.stoi_matrix is None
+
+        tlist = np.logspace(-13, -5, 81, dtype=float)
+
+        # Integrate to get the solution at each time point
+        t = []
+        y_off = []
+        species_rates_off = []
+        t.append(rxn_system.t)
+
+        # You must make a copy of y because it is overwritten by DASSL at
+        # each call to advance()
+        y_off.append(rxn_system.y.copy())
+        species_rates_off.append(rxn_system.core_species_rates.copy())
+        for t1 in tlist:
+            rxn_system.advance(t1)
+            t.append(rxn_system.t)
+            # You must make a copy of y because it is overwritten by DASSL at
+            # each call to advance()
+            y_off.append(rxn_system.y.copy())
+            species_rates_off.append(rxn_system.core_species_rates.copy())
+        run_time = time.time() - start_time
+        print(f"Simulation took {run_time:.3e} seconds")
+
+        # Convert the solution vectors to np arrays
+        t = np.array(t, float)
+        y_off = np.array(y_off, float)
+        species_rates_off = np.array(species_rates_off, float)
+
+        # Check that we've reached equilibrium
+        assert abs(species_rates_off[-1, 0] - 0.0) < 1e-2
+
+        # Check that coverages are different
+        assert not np.allclose(y, y_off)
+        assert not np.allclose(species_rates, species_rates_off)
diff --git a/test/rmgpy/thermo/convertTest.py b/test/rmgpy/thermo/convertTest.py
index 94637ce8f6..ef2f7265e3 100644
--- a/test/rmgpy/thermo/convertTest.py
+++ b/test/rmgpy/thermo/convertTest.py
@@ -57,6 +57,7 @@ def setup_class(self):
             S0=(-118.46 * constants.R, "J/(mol*K)"),
             Tmin=(10, "K"),
             Tmax=(3000, "K"),
+            thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}},
             comment="C2H6",
         )
         self.nasa = NASA(
@@ -93,6 +94,7 @@ def setup_class(self):
             E0=(-93.6077, "kJ/mol"),
             Cp0=(4.0 * constants.R, "J/(mol*K)"),
             CpInf=(21.5 * constants.R, "J/(mol*K)"),
+            thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}},
             comment="C2H6",
         )
         self.thermodata = ThermoData(
@@ -108,6 +110,7 @@ def setup_class(self):
             Tmin=(10, "K"),
             Tmax=(3000, "K"),
             E0=(-93.6077, "kJ/mol"),
+            thermo_coverage_dependence = {'OX':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}},
             comment="C2H6",
         )
 
@@ -129,6 +132,8 @@ def test_convert_wilhoit_to_nasa(self):
             s_nasa = nasa.get_entropy(T)
             assert abs(s_nasa - s_wilhoit) < 1e0
         assert abs(wilhoit.E0.value_si - nasa.E0.value_si) < 1e1
+        assert repr(nasa.thermo_coverage_dependence) != {}
+        assert repr(nasa.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence)
 
     def test_convert_wilhoit_to_thermo_data(self):
         """
@@ -149,6 +154,8 @@ def test_convert_wilhoit_to_thermo_data(self):
         s_thermodata = thermodata.get_entropy(T)
         assert round(abs(s_thermodata - s_wilhoit), 4) == 0
         assert abs(wilhoit.E0.value_si - thermodata.E0.value_si) < 1e1
+        assert repr(thermodata.thermo_coverage_dependence) != {}
+        assert repr(thermodata.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence)
 
     def test_convert_nasa_to_wilhoit(self):
         """
@@ -168,6 +175,8 @@ def test_convert_nasa_to_wilhoit(self):
             s_nasa = nasa.get_entropy(T)
             assert abs(s_nasa - s_wilhoit) < 1e0
         assert abs(nasa.E0.value_si - wilhoit.E0.value_si) < 2e1
+        assert repr(wilhoit.thermo_coverage_dependence) != {}
+        assert repr(wilhoit.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence)
 
     def test_convert_nasa_to_thermo_data(self):
         """
@@ -188,6 +197,8 @@ def test_convert_nasa_to_thermo_data(self):
         s_nasa = nasa.get_entropy(T)
         assert round(abs(s_nasa - s_thermodata), 4) == 0
         assert abs(nasa.E0.value_si - thermodata.E0.value_si) < 1e1
+        assert repr(thermodata.thermo_coverage_dependence) != {}
+        assert repr(thermodata.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence)
 
     def test_convert_thermo_data_to_wilhoit(self):
         """
@@ -208,6 +219,8 @@ def test_convert_thermo_data_to_wilhoit(self):
         s_thermodata = thermodata.get_entropy(T)
         assert round(abs(s_thermodata - s_wilhoit), 3) == 0
         assert abs(thermodata.E0.value_si - wilhoit.E0.value_si) < 1e1
+        assert repr(wilhoit.thermo_coverage_dependence) != {}
+        assert repr(wilhoit.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence)
 
     def test_convert_thermo_data_to_nasa(self):
         """
@@ -228,6 +241,8 @@ def test_convert_thermo_data_to_nasa(self):
         s_nasa = nasa.get_entropy(T)
         assert abs(s_nasa - s_thermodata) < 1e0
         assert abs(thermodata.E0.value_si - nasa.E0.value_si) < 1e1
+        assert repr(nasa.thermo_coverage_dependence) != {} 
+        assert repr(nasa.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence)
 
     def test_wilhoit_nasa_wilhoit(self):
         """
@@ -248,6 +263,8 @@ def test_wilhoit_nasa_wilhoit(self):
             s_2 = wilhoit2.get_entropy(T)
             assert abs(s_1 - s_2) < 1e0
         assert abs(wilhoit1.E0.value_si - wilhoit2.E0.value_si) < 1e1
+        assert repr(wilhoit1.thermo_coverage_dependence) != {}
+        assert repr(wilhoit1.thermo_coverage_dependence) == repr(wilhoit2.thermo_coverage_dependence)
 
     def test_wilhoit_thermo_data_wilhoit(self):
         """
@@ -268,3 +285,5 @@ def test_wilhoit_thermo_data_wilhoit(self):
             s_2 = wilhoit2.get_entropy(T)
             assert abs(s_1 - s_2) < 1e0
         assert abs(wilhoit1.E0.value_si - wilhoit2.E0.value_si) < 1e1
+        assert repr(wilhoit1.thermo_coverage_dependence) != {}
+        assert repr(wilhoit1.thermo_coverage_dependence) == repr(wilhoit2.thermo_coverage_dependence)
diff --git a/test/rmgpy/thermo/nasaTest.py b/test/rmgpy/thermo/nasaTest.py
index f08b28e0da..23255082f9 100644
--- a/test/rmgpy/thermo/nasaTest.py
+++ b/test/rmgpy/thermo/nasaTest.py
@@ -31,8 +31,7 @@
 This script contains unit tests of the :mod:`rmgpy.thermo.nasa` module.
 """
 
-import os.path
-
+import os.path, ast
 
 import numpy as np
 
@@ -69,6 +68,7 @@ def setup_class(self):
         self.Tmax = 3000.0
         self.Tint = 650.73
         self.E0 = -782292.0  # J/mol.
+        self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[(1,'J/mol'),(2,'J/mol'),(3,'J/mol')], "entropy-coefficients":[(1,'J/(mol*K)'),(2,'J/(mol*K)'),(3,'J/(mol*K)')]}}
         self.comment = "C2H6"
         self.nasa = NASA(
             polynomials=[
@@ -82,6 +82,7 @@ def setup_class(self):
             Tmin=(self.Tmin, "K"),
             Tmax=(self.Tmax, "K"),
             E0=(self.E0, "J/mol"),
+            thermo_coverage_dependence=self.thermo_coverage_dependence,
             comment=self.comment,
         )
 
@@ -136,6 +137,12 @@ def test_comment(self):
         Test that the NASA comment property was properly set.
         """
         assert self.nasa.comment == self.comment
+    
+    def test_thermo_coverage_dependence(self):
+        """
+        Test that the thermo_coverage_dependence property was properly set.
+        """
+        assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(self.thermo_coverage_dependence))
 
     def test_is_temperature_valid(self):
         """
@@ -263,6 +270,7 @@ def test_pickle(self):
         assert self.nasa.Tmax.units == nasa.Tmax.units
         assert self.nasa.E0.value == nasa.E0.value
         assert self.nasa.E0.units == nasa.E0.units
+        assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(nasa.thermo_coverage_dependence))
         assert self.nasa.comment == nasa.comment
 
     def test_repr(self):
@@ -296,6 +304,7 @@ def test_repr(self):
         assert self.nasa.Tmax.units == nasa.Tmax.units
         assert self.nasa.E0.value == nasa.E0.value
         assert self.nasa.E0.units == nasa.E0.units
+        assert ast.literal_eval(repr(self.nasa.thermo_coverage_dependence)) == ast.literal_eval(repr(nasa.thermo_coverage_dependence))
         assert self.nasa.comment == nasa.comment
 
     def test_to_cantera(self):
@@ -326,6 +335,7 @@ def test_to_nasa(self):
 
         spc = Species().from_smiles("CC")
         spc.get_thermo_data()
+        spc.thermo.thermo_coverage_dependence = self.thermo_coverage_dependence
 
         T = 1350.0  # not 298K!
 
@@ -353,6 +363,7 @@ def test_to_nasa(self):
 
         assert round(abs(s_nasa - s_td), -1) == 0
         assert wilhoit.comment == nasa.comment
+        assert repr(wilhoit.thermo_coverage_dependence) == repr(nasa.thermo_coverage_dependence)
 
         # nasa to wilhoi performed in wilhoitTest
 
@@ -364,6 +375,13 @@ def test_nasa_as_dict_full(self):
         assert nasa_dict["E0"]["value"] == self.E0
         assert nasa_dict["Tmin"]["value"] == self.Tmin
         assert nasa_dict["Tmax"]["value"] == self.Tmax
+        assert nasa_dict["thermo_coverage_dependence"].keys() == self.thermo_coverage_dependence.keys()
+        sp_name = list(self.thermo_coverage_dependence.keys())[0]
+        assert nasa_dict['thermo_coverage_dependence'][sp_name]['model'] == self.thermo_coverage_dependence[sp_name]['model']
+        enthalpy_list = nasa_dict['thermo_coverage_dependence'][sp_name]['enthalpy-coefficients']['object'] 
+        # return [(str(coeff.value), str(coeff.units))for coeff in enthalpy_list], self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients']
+        assert [(int(coeff.value), str(coeff.units))for coeff in enthalpy_list] == self.thermo_coverage_dependence[sp_name]['enthalpy-coefficients']
+        assert nasa_dict['thermo_coverage_dependence'][sp_name]['entropy-coefficients']['object'] == self.thermo_coverage_dependence[sp_name]['entropy-coefficients']
         assert nasa_dict["comment"] == self.comment
         assert tuple(nasa_dict["polynomials"]["polynomial1"]["coeffs"]["object"]) == tuple(self.coeffs_low)
         assert tuple(nasa_dict["polynomials"]["polynomial2"]["coeffs"]["object"]) == tuple(self.coeffs_high)
diff --git a/test/rmgpy/thermo/thermodataTest.py b/test/rmgpy/thermo/thermodataTest.py
index e3a633c2da..3c6b39e7e0 100644
--- a/test/rmgpy/thermo/thermodataTest.py
+++ b/test/rmgpy/thermo/thermodataTest.py
@@ -53,6 +53,7 @@ def setup_class(self):
         self.Tmin = 100.0
         self.Tmax = 3000.0
         self.E0 = -782292.0
+        self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}
         self.comment = "C2H6"
         self.thermodata = ThermoData(
             Tdata=(self.Tdata, "K"),
@@ -64,6 +65,7 @@ def setup_class(self):
             Tmin=(self.Tmin, "K"),
             Tmax=(self.Tmax, "K"),
             E0=(self.E0, "J/mol"),
+            thermo_coverage_dependence=self.thermo_coverage_dependence,
             comment=self.comment,
         )
 
@@ -130,6 +132,12 @@ def test_comment(self):
         Test that the ThermoData comment property was properly set.
         """
         assert self.thermodata.comment == self.comment
+    
+    def test_thermo_coverage_dependence(self):
+        """
+        Test that the ThermoData thermo_coverage_dependence property was properly set.
+        """
+        assert repr(self.thermodata.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence)
 
     def test_is_temperature_valid(self):
         """
@@ -261,6 +269,7 @@ def test_pickle(self):
         assert round(abs(self.thermodata.E0.value - thermodata.E0.value), 4) == 0
         assert self.thermodata.E0.units == thermodata.E0.units
         assert self.thermodata.label == thermodata.label
+        assert repr(self.thermodata.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence)
         assert self.thermodata.comment == thermodata.comment
 
     def test_repr(self):
@@ -295,6 +304,7 @@ def test_repr(self):
         assert round(abs(self.thermodata.E0.value - thermodata.E0.value), 4) == 0
         assert self.thermodata.E0.units == thermodata.E0.units
         assert self.thermodata.label == thermodata.label
+        assert repr(self.thermodata.thermo_coverage_dependence) == repr(thermodata.thermo_coverage_dependence)
         assert self.thermodata.comment == thermodata.comment
 
     def test_is_all_zeros(self):
diff --git a/test/rmgpy/thermo/wilhoitTest.py b/test/rmgpy/thermo/wilhoitTest.py
index 4b74fd72a7..6a350d7833 100644
--- a/test/rmgpy/thermo/wilhoitTest.py
+++ b/test/rmgpy/thermo/wilhoitTest.py
@@ -45,7 +45,6 @@ class TestWilhoit:
     """
     Contains unit tests of the :class:`Wilhoit` class.
     """
-
     def setup_class(self):
         self.Cp0 = 4.0
         self.CpInf = 21.5
@@ -59,6 +58,7 @@ def setup_class(self):
         self.Tmin = 300.0
         self.Tmax = 3000.0
         self.comment = "C2H6"
+        self.thermo_coverage_dependence = {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}':{'model':'polynomial', 'enthalpy-coefficients':[1,2,3], "entropy-coefficients":[1,2,3]}}
         self.wilhoit = Wilhoit(
             Cp0=(self.Cp0 * constants.R, "J/(mol*K)"),
             CpInf=(self.CpInf * constants.R, "J/(mol*K)"),
@@ -71,6 +71,7 @@ def setup_class(self):
             S0=(self.S0 * constants.R, "J/(mol*K)"),
             Tmin=(self.Tmin, "K"),
             Tmax=(self.Tmax, "K"),
+            thermo_coverage_dependence=self.thermo_coverage_dependence,
             comment=self.comment,
         )
 
@@ -159,6 +160,12 @@ def test_comment(self):
         Test that the Wilhoit comment property was properly set.
         """
         assert self.wilhoit.comment == self.comment
+    
+    def test_thermo_coverage_dependence(self):
+        """
+        Test that the Wilhoit thermo_coverage_dependence property was properly set.
+        """
+        assert repr(self.wilhoit.thermo_coverage_dependence) == repr(self.thermo_coverage_dependence)
 
     def test_is_temperature_valid(self):
         """
@@ -287,6 +294,7 @@ def test_pickle(self):
         assert self.wilhoit.Tmax.units == wilhoit.Tmax.units
         assert round(abs(self.wilhoit.E0.value - wilhoit.E0.value), 4) == 0
         assert self.wilhoit.E0.units == wilhoit.E0.units
+        assert repr(self.wilhoit.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence)
         assert self.wilhoit.comment == wilhoit.comment
 
     def test_repr(self):
@@ -318,6 +326,7 @@ def test_repr(self):
         assert self.wilhoit.Tmax.units == wilhoit.Tmax.units
         assert round(abs(self.wilhoit.E0.value - wilhoit.E0.value), 1) == 0
         assert self.wilhoit.E0.units == wilhoit.E0.units
+        assert repr(self.wilhoit.thermo_coverage_dependence) == repr(wilhoit.thermo_coverage_dependence)
         assert self.wilhoit.comment == wilhoit.comment
 
     def test_fit_to_data(self):
@@ -381,6 +390,7 @@ def test_to_wilhoit(self):
 
         spc = Species().from_smiles("CC")
         spc.get_thermo_data()
+        spc.thermo.thermo_coverage_dependence = self.thermo_coverage_dependence
 
         T = 1350.0  # not 298K!
 
@@ -448,6 +458,11 @@ def test_wilhoit_as_dict(self):
                 "value": 178.76114800000002,
             },
             "class": "Wilhoit",
+            'thermo_coverage_dependence': {'1 O u0 p2 c0 {2,D} \n 2 X u0 p0 c0 {1,D}': {
+                                                                                        'model': 'polynomial', 
+                                                                                        'enthalpy-coefficients': {'class': 'np_array', 'object': [1, 2, 3]}, 
+                                                                                        'entropy-coefficients': {'class': 'np_array', 'object': [1, 2, 3]}}
+                                          }
         }
 
     def test_make_wilhoit(self):
@@ -456,6 +471,6 @@ def test_make_wilhoit(self):
         """
         wilhoit_dict = self.wilhoit.as_dict()
         new_wilhoit = Wilhoit.__new__(Wilhoit)
-        class_dictionary = {"ScalarQuantity": ScalarQuantity, "Wilhoit": Wilhoit}
+        class_dictionary = {"ScalarQuantity": ScalarQuantity, "Wilhoit": Wilhoit, "np_array": np.array}
 
         new_wilhoit.make_object(wilhoit_dict, class_dictionary)