We present an interactive 3-D computer program (LitMod3D) developed to perform combined geophysical-petrological modeling of the lithosphere and sublithospheric upper mantle. In contrast to other available modeling software, LitMod3D is built within an internally consistent thermodynamicgeophysical framework, where all relevant properties are functions of temperature, pressure, and composition. By simultaneously solving the heat transfer, thermodynamic, rheological, geopotential, and isostasy (local and flexural) equations, the program outputs temperature, pressure, surface heat flow, density (bulk and single phase), seismic wave velocities, geoid and gravity anomalies, elevation, and lithospheric strength for any given model. These outputs can be used to obtain thermal and compositional models of the lithosphere and sublithospheric upper mantle that simultaneously fit all available geophysical and petrological observables. We illustrate some of the advantages and limitations of LitMod3D using synthetic models and comparing our predictions with those from other modeling methods. In particular, we show that (1) temperature at midlithosphere depths may be overestimated by as much as 200 K when compositional heterogeneities in the mantle and T-P effects are not considered in lithospheric models and (2) the neglect of mantle phase transformations on gravity-based models in thin-crust settings can result in a significant overestimation and underestimation of the derived crustal thickness and its internal density distribution, respectively.