In the present paper, mathematical modeling and numerical solution for the 3D carbonation phenomenon on cement mortars involving aggregates have been investigated for the first time. To achieve this assessment, the Papadakis analytical proposal has been fully investigated under 1D and 3D considerations. The molar concentration variations of the hydrates (CSH and Ca(OH) 2 ) and unhydrated grains (C 2 S and C 3 S) have been analyzed during carbonation, including hydration phenomenon. At the first stage, 1D numerical modeling of multi-reaction carbonation has been applied to atmospheric carbonation and verified by experiments in the literature. Afterwards, numerical simulations have been achieved on 3D numerical mortar samples for accelerated carbonation including aggregates using the relevant granulometry, whose applications sustain more realistic outcomes. The numerical solution has been done using the finite element method for the extremely non-linear transient system of PDEs. The numerical experiments have been compared to those done using the pH detector (Phenolphthalein and Thymophtaleine) and Differential Thermal Analysis (DTA). Some conclusions and outlooks pertaining to carbonation modeling have been emphasized, including water migration and CaCO 3 precipitation issues in improving the analytical modeling for accelerated carbonation cases.