A robust computational approach is proposed and validated for reconstructing constitutive relations (material properties) in complex multiphysics phenomena based on incomplete and noisy measurements. The method is applicable to broad range of problems arising in nonequilibrium thermodynamics and continuum mechanics where, in contrast to “classical” parameter estimation problems, the control variable is a function of state (i.e. the dependent variable) rather than the independent variable. An elegant and computationally efficient solution for this inverse problem is obtained by formulating it as a PDE-constrained optimization problem which can be solved with a gradient-based technique within the “optimize-then-discretize” framework. The main novelty is that the nonlinear constitutive relation is determined in a very general form with no prior knowledge except some assumptions on smoothness. A key element of the proposed approach is the cost functional gradients are expressed in terms of integrals defined over the level sets of the state variable field. Robustness of the method is shown in applications for two different classes of problems, namely identification of temperature-dependent material properties in complex thermo-fluid phenomena and identification of inertial manifolds in reduced-order models of hydrodynamic instabilities.