Spatial discretization is critical for robustly modeling spatially distributed hydrologic processes, particularly runoff routing. Flow routing models that use a uniform-resolution cartesian grid have several limitations including inconstancy of travel time in different directions, inaccurate representation of watershed boundary with sharp corners, and lacking of interface connection between land and ocean in Earth Systems Models (ESMs). The different types of grids used by the river, land and ocean components in ESMs lead to significant challenges in capturing the river-land-ocean continuum. Earlier studies have suggested that the use of a hexagonal grid within flow routing models has the potential to resolve the aforementioned limitations, but applications of such grids are rare in ESMs. In this study, we extend MOSART, the flow routing model of the Energy Exascale Earth System Model (E3SM), to use a hexagonal grid. We evaluate MOSART simulations that use hexagonal and cartesian grids against multiple observational datasets and compare the performance at multiple spatial resolutions. This study improves our understanding of the impacts of spatial discretization on flow routing model performance and the corresponding uncertainties. It also paves the way to better coupling river, land, and ocean components in ESMs.