Turbulence closure models will continue to be necessary in order to perform computationally affordable simulations in the foreseeable future. It is expected that Reynolds-averaged Navier-Stokes (RANS) turbulence models will still be useful with the further development of the more accurate, but computationally expensive large eddy simulation (LES), especially in industry. The use of the robust but often inaccurate linear eddy viscosity closures is still widespread in industry. More complex closure models, such as Reynolds stress models and nonlinear eddy viscosity models, provide a more general description of the underlying physics of turbulent flows. Nevertheless, because of implementational difficulties or failure to provide consistent improvements over the more robust linear models, RANS turbulence modeling is considered to have reached a plateau. In the past few years, the availability of high-fidelity datasets, the increased accuracy of machine learning algorithms, and the rise in computational power led to the proposal of several data-driven approaches to turbulence modeling. The general idea is to use experimental and high-fidelity data to develop or enhance RANS turbulence models, instead of employing an approach purely based on physics. As in any emerging field, there are many possibilities for further developing the novel approaches to data-driven turbulence modeling. Recent work combined machine learning with statistical inversion. First, a spatially varying correction is applied to the RANS model and optimized by minimizing the discrepancy between the RANS output and the data for several flows. Machine learning is used to approximate a function between a set of flow features and the inferred corrections. The aim of this work is to further investigate this methodology, called the paradigm of field inversion and machine learning, in a broader set of test cases by inferring a spatially varying correction to the production term of the ω-equation in the k−ω model and to the eigenvalues of the Reynolds stress tensor. The gradients of this high-dimensional optimization process are obtained by implementing the continuous adjoint of the k−ω model in OpenFOAM. Gaussian processes and random forests are then used to approximate a function between mean flow features and the inferred corrections. It was found that both formulations are able to accurately infer the mean velocities and related quantities of interest, but that the inferred corrective terms are often non-unique or not physically interpretable. For several flow cases, the corrective terms were able to generalize to unseen Reynolds number and flow geometries.