Acidizing is a widely used technique for intensification of geothermal wells. At certain conditions, acid injection results in a formation of highly conductive channels, e.g. wormholes. The associated rock matrix dissolution is highly nonlinear and may influence the fluid flow as well as rock properties. Wormhole patterns and breakthrough time are mainly controlled by acid strength and injection rate. Due to high complexity and, therefore, high computational cost, existing conventional simulators usually do not include complex chemical reactions or use weak (sequential) coupling of flow with reactive transport. Although the sequential approach is good for various applications, it is limited by the CFL condition and might suffer from poor convergence. In this project, we investigate the wormhole phenomenon to provide a stable fully implicit method (FIM) to numerically solve the reactive flow and transport problem coupled with equilibrium and kinetic reactions associated with the acid injection. The developed framework is aimed to predict a pattern of dissolution as well as wormhole breakthrough time. For an accurate description of species interactions, we directly connect PHREEQC equilibrium computation with FIM framework for kinetic dissolution using the Element Balance approach. This coupling was achieved with an application of a recently introduced Operator-Based Linearization scheme utilized within a Delft Advanced Research Terra Simulator (DARTS).