Knowledge about the Earth's stress field and its sources can provide better understanding and interpretation of geodynamic and tectonic processes and regimes in the Earth's lithosphere. Stresses can be measured with different in-situ techniques and analysed by the study of focal mechanisms and stress induced geological structures. Quantifying single stress sources however remains a difficult and not uncommonly vague procedure. Modelling stress contributions can provide principle insight into the total and relative stress composition. One contribution is the internal gravitational stress in the lithosphere, induced by lateral density variation. The leading quantity is the Geopotential Energy, the integrated lithostatic pressure in a rock column, which is related to horizontal stresses by the Equations of Equilibrium. The Geopotential Energy can be furthermore linearly related to the Geoid under assumption of local isostasy. Satellite Geoid measurements contain, however, deeper mantle responses of most likely longwavelength. Still after filtering, the Geoid can't be satisfyingly corrected. Existing shallow signals can be hereby extinguished as well, for instance the somewhat age dependent signal of the oceanic lithosphere. An entire modelling of the shallow Geopotential Energy is hereby approached, not taking into account possible deeper signals but all lithospheric signals for the subsequent stress calculation. Therefore a global lithospheric density model is necessary to calculate the corresponding response to Geopotential Energy and the Geoid. A linearized inverse method fits a lithospheric reference model to reproduce measured data sets, such as topography and surface heat flow, while assuming isostasy and solving the steady state heat equation. A FEM code solves the equations of equilibrium of stresses for a three dimensional spherical shell. The modelled results are shown and compared with concurrent publications.