Reactive transport problems involve the coupling between the chemical interactions of different species and their transport by advection and diffusion. It leads to the solution of a non-linear systems of partial differential equations coupled to local algebraic or differential equations. Developing software for these two components involves fairly different techniques, so that methods based on loosely coupled modules are desirable. On the other hand, numerical issues such as robustness and convergence require closer couplings, such as simultaneous solution of the overall system. The method described in this paper allows a separation of transport and chemistry at the software level, while keeping a tight numerical coupling between both subsystems. We give a formulation that eliminates the local chemical concentrations and keeps the total concentrations as unknowns, then recall how each individual subsystem can be solved. The coupled system is solved by a Newton–Krylov method. The block structure of the model is exploited both at the nonlinear level, by eliminating some unknowns, and at the linear level by using block Gauss-Seidel or block Jacobi preconditioning. The methods are applied to a 1D case of the MoMaS benchmark.