Modeling　reactive　transport　in　discontinuous　and　heterogeneous　porous　media　is　key　to　the　understanding　of　geochemical　systems.　Integral　formulations　of　conservation　equations　can　be　an　alternative　approach　to　solving　transport　problems　in　such　complex　media　compared　to　the　classical　local　(differential)　formulations.　A　prominent　example　is　Peridynamics.　However,　it　has　been　developed　only　for　advection-dispersion　transport　with　a　simple　bimolecular　reaction,　which　is　of　limited　utility　in　complex　geochemical　problems　in　real　environments.　This　study　integrates　bond-based　Peridynamics　with　advanced　geochemical　modeling　(PHREEQC),　enabling　accurate　simulation　of　reactive　transport　in　heterogeneous　porous　systems,　overcoming　limitations　in　grid-based　discretization　methods.　The　sequential　non-iterative　approach　is　introduced　to　address　the　coupling　between　the　transport　of　chemical　species　(solved　by　Peridynamics)　and　geochemical　reactions　(solved　by　PHREEQC).　The　proposed　model　is　verified　with　a　set　of　benchmarks.　A　series　of　cases　is　studied　to　show　the　model＇s　capabilities　in　the　prediction　of　reactive　transport　in　porous　media　with　fractures　and　heterogeneities　(e.g.,　permeable　and　impermeable　inclusions).　The　current　model　can　be　used　to　quantify　the　impact　of　the　permeable/non-permeable　inclusions　on　non-uniform　solution　migration　and　sharp　fronts　of　mineral　precipitation/dissolution　without　any　refinements　and　modifications　at　the　interface.　The　PD　reactive　transport　also　captures　the　influences　of　fractures　in　accelerating　the　solute　transport　and　enhancing　the　mineral　reaction　rate.　The　current　multi-physics　nonlocal　reactive　transport　formulations　can　be　easily　extended　to　study　more　complex　problems,　such　as　reactive　flow-mechanical　process　coupled　behavior　and　accurate　description　of　the　mineral　dissolution/precipitation　interface　at　the　micro-scale.