An automated method is presented to identify energetically feasible mechanisms of heterogeneous catalytic systems comprising of several thousand reactions and species. Specifically, it combines automated rule-based network generation of complex networks with semi-empirical estimation of thermochemical properties to (a) generate large reaction networks, (b) calculate thermochemistry and activation barriers of each step on-the-fly, and (c) identify energetically feasible pathways and mechanisms to experimentally observed products. This method is applied to analyze the mechanism of glycerol conversion on transition metal catalysts for which a network of 3300 reactions and 500 species was generated. Using (i) group additivity methods for calculating reaction and species energies, (ii) linear free energy relationships for activation barriers, and (iii) linear scaling relationships for comparative assessment of energetics of species on different transition metals, the energetically feasible pathways for forming C-C scission and C-O scission products (syngas and propane diol, respectively) are identified on platinum, palladium, rhodium, and ruthenium catalysts. Our results indicate that (a) C-C scission products are preferred on platinum and palladium, (b) rhodium and ruthenium will have a comparatively higher selectivity for C-O scission products, and (c) glycerol tends to undergo several dehydrogenation steps prior to undergoing C-C scission, which are all consistent with experimental observations and DFT calculations. We propose that this method can be used to screen a large number of pathways in complex catalytic networks and to thereby identify the dominant modes in the system, because it is fast and scalable in handling larger networks, flexible in accommodating different types of semi-empirical correlations, and generic in handling any catalytic chemistry.