<p>Superlinear convergence occurs frequently as a desirable second stage in a Krylov iteration, and its understanding has undergone a deep development. This paper aims to contribute to this long history of analysis for a class of nonlinear systems of partial differential equations (PDEs) covering reaction–diffusion–convection processes, discretized with finite elements (FE), linearized with a Newton–Krylov method and applying suitable operator preconditioning. We give practical estimates of the rate of superlinear convergence of the arising Krylov iterations in terms of the accessible data of the PDEs. We are in particular interested in the effect of integrability properties of possibly unbounded sources. We obtain robust superlinear rates both in the sense of mesh independence and of uniform behaviour w.r.t. the outer Newton iterations. The numerical examples reinforce the theoretical results.</p>