In this study, an effective method is developed to couple the multi-directional pore-network (MDPN) model, called PoreFlow, with the chemistry simulator PHREEQC in order to simulate the reactive flow of chelating agents into carbonate media. A novel method based on solving a 1D chemical reactions and projecting the chemistry in the 3D pore scale flow field is developed, which significantly increased the computational efficiency. Results from several simulations are presented to demonstrate the capability of PoreFlow and the developed method to mimic reactive transport of chelating agents in porous media. It is also shown that the interactions between chelating agents and carbonate rocks can be accurately modeled by implementing the proper multi-step protonation kinetic of chelating agents versus pH. The developed model is capable of computing the change in porosity and the 3D flow field due to calcite dissolution by the progress of chemical reactions. In addition, the evolving flow model provides the corresponding relationship between porosity and permeability for a specific porous media. The model successfully reproduced the experimentally observed permeability changes due to the creation of wormholes at various acid flow rates. The results are validated by comparing the simulation outputs to the micro-CT images taken from the samples after the chelating experiments. Moreover, it is shown that the average dissolution rate for heterogeneous carbonate rocks is not only governed by the pore-scale physical heterogeneities but also by the local chemical equilibrium and speciation of chelating agents at various pH values.