Hydrodynamic models are essential for understanding water bodies’ circulation, and they complement studies regarding pollutant transport, especially concerning estuarine and river systems, e.g., the ones surrounding the city of Belém in northern Brazil. Hydrodynamic simulations were developed, using Guamá river (GR), Guajará bay (GB), and Pará river estuary (PRE), to assess circulation patterns under high (HD) and low river discharge (LD) and along tidal cycles. The simulations were run in Delft3D-FM, using an unstructured grid with resolutions from 1,000,000 to 40 m². Bathymetry was interpolated from observed data, while Manning’s roughness ranged from 0.020 to 0.039. Tidal boundaries were provided by the TPXO 7.2 Inverse Tidel Model—PRE mouth—and in situ observations—Breves Strait—while discharge boundaries were applied to the Guamá, Acará, Tocantins and Moju rivers, supplied with values estimated from historical records of the Brazilian National Waters Agency. Model validation, based on Pearson’s coefficient, NSE, and pRMSE, indicated efficient performance. Under HD, PRE showed higher velocities on the lower estuary during flood tides and on the higher estuary during the ebb, with persistent upstream flux driven by the Tocantins river discharge. In GB and GR, velocities ranged between 0.5 and 1.0 m/s during flood tide, with a low-dynamics zone in GR, near Belém, due to opposing flows from GB and GR, which disappeared throughout the ebb. Under LD, PRE’s patterns remained similar, though velocities slightly increased on the higher estuary during flood and on the lower estuary during ebb. In GB, velocities peaked during flood along the left margin due to a deep channel, while low-dynamic zones disappeared. Overall, PRE exhibited stable circulation whilst GR and GB differed between season and tidal cycles. These findings can highlight critical zones for pollutant retention and dispersion, providing a basis for water quality assessments in the region.