Abstract: Dissolution-driven fracture governs the long-term integrity of geological CO2 storage, geothermal reservoirs, and nuclear waste repositories. Modeling it requires solving strongly coupled mechanical, chemical, and phase-field equations. The solver choice controls their accuracy, cost, and convergence. This study investigates the role of solution strategies in dissolution-driven fracture in geomaterials. The coupled model is implemented in Abaqus user-element (UEL) subroutines, with all tangent blocks derived analytically. Five solvers are compared: (1) block-diagonal Newton monolithic, (2) BFGS quasi-Newton monolithic, (3) single-pass fully staggered, (4) single-pass partially staggered, and (5) hybrid multi-pass fully staggered. Performance is quantified by accuracy (errors in peak load, peak and failure displacements, and local CO2 concentration), cost (CPU and wall-clock times), and robustness, measured by new factors that assess convergence over the whole simulation and capture iteration spikes. The solvers are tested on mechanical and chemo-mechanical tension, shear, and multiple-crack benchmarks, including CO2 diffusivities raised up to a thousandfold. The results reveal a clear trade-off between accuracy, cost, and robustness. The quasi-Newton solver matches the Newton reference at a fraction of its cost but can fail under intense transport. The single-pass fully staggered solver is inexpensive and always completes, but its post-peak response depends strongly on the time step. The hybrid multi-pass solver also always completes, retains near-reference accuracy, and avoids iteration spikes at moderate cost. The partially staggered solver is the least robust. These findings guide practical selection: the single-pass staggered solver suits rapid screening, quasi-Newton offers accuracy when convergent, and hybrid multi-pass provides the best balance.