Reservoir stimulation is a widely used technique in the oil and gas industry for increasing the productivity of hydrocarbon reservoirs, most notably in carbonate formations. This work aims to develop an optimization workflow under uncertainty for matrix acidizing. A reactive transport model is implemented in a finite-element framework to simulate the initiation and propagation of dissolution channels in porous carbonate rock. The model is verified using an analytical solution. We utilize surrogate modeling based on polynomial chaos expansion (PCE) and Sobol indices to identify the most significant parameters. We investigate the effect of varying 12 identified parameters on the efficiency of the stimulation process using dimensionless groups, including the Damköhler, Peclet, and acid capacity numbers. Furthermore, the surrogate model reproduces the physics-based results accurately, including the dissolution channels, the pore volume to breakthrough, and the effective permeability of the stimulated rock. The developed workflow assesses how uncertainties propagate to the model's response, where the surrogate model is used to calculate the univariate effect. The global sensitivity analysis shows that the acid capacity number is the most significant parameter for the pore volume to breakthrough with the highest Sobol index. The marginal effect calculated for the individual parameter confirms the results from Sobol indices. This work provides a systematic workflow for uncertainty analysis and optimization applied to the processes of rock stimulation. Characterizing the impact of uncertainty provides physical insights and a better understanding of the matrix acidizing process.