In numerical codes for reactive transport modeling, systems of nonlinear chemical equations are often solved through the Newton Raphson method (NR). NR is an iterative procedure that results in a sequential solution of linear systems. The algorithm is known for its effectiveness in the vicinity of the solution but also for its lack of robustness otherwise. Therefore, inaccurate initial conditions can lead to non-convergence or excessive numbers of iterations, which significantly increase the computational cost. In this work, we show that inaccurate initial conditions can lead to very ill-conditioned system matrices, which makes NR inefficient. This efficiency is improved by preconditioning techniques and/or by coupling the NR method with a zero-order method called the positive continuous fraction method. Numerical experiments that are based on seven different test cases show that the ill-conditioned linear systems within NR represent a problem and that coupling NR with a method that bypasses the computation of the Jacobian matrix significantly improves the robustness and efficiency of the algorithm. VC 2016 American Institute of Chemical Engineers AIChE J, 63: 1246-1262, 2017