Accurate surrogate modeling in engineering is often constrained by the high computational cost of generating training data from large-scale numerical simulations. In many industrial applications, only a limited number of simulations can be afforded, which severely restricts the achievable surrogate accuracy, particularly in high dimensional parametric spaces. Advances in differentiable simulation and adjoint methods enable the efficient computation of first- and higher order sensitivities of simulation outputs with respect to input parameters. This development opens new opportunities for Sobolev training of neural networks, where derivative information provides valuable local structure that can significantly enhance data efficiency. However, Sobolev training is inherently challenging in practice, as it requires reliable and scalable access to derivative information, especially for complex numerical solvers. While adjoint-based derivatives have been successfully employed in first order gradient enhanced kriging, their systematic use for training neural network surrogates of finite element (FE) and computational fluid dynamics (CFD) simulations remains largely unexplored. Existing approaches to second-order Sobolev training in finite element settings rely on analytically derived gradients and Hessians, which substantially limits their applicability to realistic industrial simulations with complex governing equations and discretizations. In this work, we propose a universal framework for Sobolev based surrogate modeling that leverages adjoint calculations from numerical simulations. Our approach is solver-agnostic for differentiable FE and CFD models, enabling gradient- and Hessian enhanced learning under strict computational budgets. The resulting framework is particularly well suited for high dimensional parametric problems, where derivative information provides a powerful mechanism for improving data efficiency. The proposed methods are evaluated on the Rosenbrock function as well as on representative industrial examples using a differentiable finite element solver. We investigate neural network training with both first- and second order Sobolev formulations, demonstrating the benefits of adjoint-based derivative information across analytical benchmarks and realistic simulation driven scenarios.