Hydrogen (H2) production has garnered attention amongst renewable energy researchers and intergovernmental agencies because of its potential to be an energy carrier and ability to satisfy the energy requirements of society. To date, the primary resource for H2 production is methane (CH4) found in natural gas, and to reduce fossil fuel utilization, other hydrocarbon sources have been identified. Glycerol is the by-product of biodiesel production, and with excess crude glycerol on the market, it proves advantageous to apply reforming methods to this feedstock to produce hydrogen. This work presents numerical models for a catalytic auto-thermal reforming (ATR) process developed for converting synthetic crude glycerol (CG) to hydrogen gas in packed bed tubular reactor (PBTR). The models were heterogenous, considering the presence of the solid catalyst and created with two numerical methods, the finite difference method (FDM) in MATLAB and the finite element method (FEM) in COMSOL Multiphysics. Both models were validated with experimental data obtained from a laboratory-scale ATR process in the presence of a Ni-based catalyst, developed under the Advanced Green Energy Systems research portfolio at the University of Regina. The power law rate model derived from the experiments was used in the model development. The CG conversion in the simulation data from the models were in good agreement with experimental data giving absolute average deviation(AAD) values of 7.56% and 6.34% for that created with FEM and FDM respectively. The conversion and temperature trends for the two models were compared to each other, with a previously built pseudo- homogeneous model for this same ATR process, developed in other work. From a qualitative standpoint, the two heterogeneous models, though not aligned perfectly, trended in the same manner. At the reaction zone, the FDM model prediction gave a slightly higher bulk fluid temperature. The heterogeneous models predicted a lower reaction zone temperature by approximately 6 oC, thus showing the impact of inclusion of the heat and mass transfer at the fluid-solid interface between into the governing equations. A different approach was used in COMSOL by manipulating the modules to overcome the limitation of generating temperature profiles for the catalyst pellet. The generated profiles for the catalyst particles in FDM and FEM were also showed to be in good agreement. Simulations were done with the FEM model varying the feed temperature from 773K to 923K and the space time from 12.71 g cat/min mol C to 158.23 g cat/min mol C. An increase in temperature resulted in increased CG conversion and H2 yield. This was also observed as the space time increased. Effectiveness factors were calculated using concentration and temperature profiles for the solid catalyst generated with the FDM to quantify how much of the catalyst surface is being used for the reaction. Initial simulations were done in accordance with the physical experiments at a catalyst size of 0.8 mm. Increasing the catalyst size in the simulation showed a significant difference between bulk fluid and solid catalyst temperature and concentration as species travelled along the reactor, decreasing CG conversion at the exit of the reactor, highlighting the importance of reducing the resistance to heat and mass resistance from fluid to solid.