Lateral density and topography variations yield in and important contribution to the lithospheric stress field. The leading quantity is the Geopotential Energy, the integrated lithostatic pressure in a rock column. The horizontal gradient of this quantity is related to horizontal stresses through the Equations of equilibrium of stresses. The Geopotential Energy furthermore can be linearly related to the Geoid under assumption of local isostasy. Satellite Geoid measurements contain, however, also non-isostatic deeper mantle responses of long wavelength. Unfortunately, high-pass filtering of the Geoid does not suppress only the deeper sources. The age-dependent signal of the oceanic lithosphere, for instance, is of long wave length and a prominent representative of in-plane stress, derived from the horizontal gradient of isostatic Geoid anomalies and responsible for the ridge push effect. Therefore a global lithospheric density model is required in order to isolate the shallow Geoid signal and calculate the stress pattern from isostatically compensated lithospheric sources. We use a linearized inverse method to fit a lithospheric reference model to observations such as topography and surface heat flow in the presence of local isostasy and a steady state geotherm. Subsequently we use a FEM code to solve the Equations of equilibrium of stresses for a three dimensional elastic shell. The modelled results are shown and compared with the global stress field and other publications.