River routing is a critical component of Earth system models (ESMs). Most existing river routing models within ESMs utilize a Cartesian rectangle mesh system, which has several limitations. For instance, it cannot be seamlessly integrated with the unstructured mesh-based ocean component. Additionally, regular discretization fails to capture drainage connectivity effectively because it does not adapt the mesh to river networks and other hydrologic features. In this study, we extend MOSART, the flow routing model of the DOE's Energy Exascale Earth System Model (E3SM), to an unstructured mesh. We apply the model to river basins at different latitudes using kilometer-scale unstructured Voronoi MPAS meshes. We evaluate MOSART's performance using both Cartesian rectangle-based datasets and observational datasets. This study enhances our understanding of the impacts of spatial discretization, particularly at the kilometer scale, on flow routing model performance and the associated uncertainties. It also lays the groundwork for improved coupling of land, river, lake, and ocean components in ESMs.