In the last few decades, hydrogeochemical problems have benefited from the strong interest in numerical modeling. One of the most recognized hydrogeochemical problems is the dissolution of the calcite in the mixing zone below limestone coastal aquifer. In many works, this problem has been modeled using a coupling algorithm between a density-dependent flow model and a geochemical model. A related difficulty is that, because of the high nonlinearity of the coupled set of equations, high computational effort is needed. During calcite dissolution, an increase in permeability can be identified, which can induce an increase in the penetration of the seawater into the aquifer. The majority of the previous studies used a fully coupled reactive transport model in order to model such problem. Romanov and Dreybrodt (J Hydrol 329:661-673, 2006) have used an alternative approach to quantify the porosity evolution in mixing zone below coastal carbonate aquifer at steady state. This approach is based on the analytic solution presented by Phillips (1991) in his book Flow and Reactions in Permeable Rock, which shows that it is possible to decouple the complex set of equation. This equation is proportional to the square of the salinity gradient, which can be calculated using a density driven flow code and to the reaction rate that can be calculated using a geochemical code. In this work, this equation is used in nonstationary step-by-step regime. At each time step, the quantity of the dissolved calcite is quantified, the change of porosity is calculated, and the permeability is updated. The reaction rate, which is the second derivate of the calcium equilibrium concentration in the equation, is calculated using the PHREEQC code (Parkhurst and Apello 1999). This result is used in GEODENS (Bouhlila 1999; Bouhlila and Laabidi 2008) to calculate change of the porosity after calculating the salinity gradient. For the next time step, the same protocol is used but using the updated porosity and permeability distributions.