The possible ill conditioning of the Reynolds average Navier–Stokes (RANS) equations when an explicit data-driven Reynolds stress tensor closure is employed is a discussion of paramount importance. This matter has far-reaching consequences on the emerging field of data-driven turbulence modelling, as well as in Reynolds stress models and in epistemic uncertainty quantification. In the present work, we explore fundamental aspects of this problem with the aid of direct numerical simulation (DNS) databases of the turbulent flows in a plane channel, in a square duct and in periodic hill geometries. We show that the RANS equations are ill conditioned in the whole range of cases analysed, even when the linear term of the Reynolds stress tensor is treated implicitly, when no information with respect to the DNS mean velocity field is provided. That is, in more than ten different simulations varying Reynolds number and geometry, using the DNS Reynolds stress tensor solely, which carries a small error, the error propagation to the mean velocity field is significantly amplified. Accurate solutions can be obtained with implicit or explicit procedures that include information from the DNS mean velocity field. We propose a new strategy along these lines for solving the RANS equations that combines ideas put forward in the literature and show that this new procedure outperforms the ones previously adopted to mitigate the error propagation to the mean velocity field. We have shown advantages of adopting the Reynolds force vector, the divergence of the Reynolds stress tensor, as a quantity with its own identity, which can impact on the choice of the target quantity to be modelled in RANS equations.