In Earth System Models, river models traditionally use structured grids (e.g., 0.5° × 0.5°) for numerical simulations. While these grids are widely applied, they have the following limitations: (1) the use of fixed resolutions cannot accurately represent real-world river networks; (2) the ability to simulate interactions between the river and ocean at the coastline is hindered as ocean models often utilize unstructured meshes; (3) the grid deformation intensifies at higher latitudes, resulting in greater model uncertainty in high-latitude regions. We propose a novel river modeling approach based on variable-resolution unstructured meshes to address these challenges. This method is primarily developed using the unstructured Voronoi mesh used in the Model for Prediction Across Scales (MPAS) of the Energy Exascale Earth System Model (E3SM) and a new flow direction model. The approach offers several advantages, including (1) seamless connectivity between land, ocean, and rivers; (2) improved spatial resolution ranging from 2 to 10km, with the entire coastline resolved; and (3) high-fidelity representation of rivers and dam locations. We apply this method to tropical, temperate, and Arctic river basins that are both rainfall and snowmelt dominated for numerical simulations and compare the results with traditional approaches. The results demonstrate that the new method achieves comparable performance to traditional methods while enhancing spatial resolution. This method provides a new approach and foundation for further improvements in land surface hydrology for Earth system modeling.