To model glacial triggering of earthquakes, it is necessary to obtain the spatio-temporal variation of glacial isostatic adjustment-induced stress during a glacial cycled. This can be computed efficiently using commercial Finite Element codes with appropriate modifications to include the important effects of ‘pre-stress advection’, ‘internal buoyancy’ and ‘self-gravity’. The modifications described in Wu (2004) are reviewed for incompressible and so-called materially compressible flat-earths. When the glacial isostatic adjustment-induced stress is superimposed on the background tectonic stress and overburden pressure, the time variation of earthquake potential at various locations in the Earth can be evaluated for any fault orientation. To model more complex slip and fault behavior over time, the three-stage Finite Element model approach of Steffen et al. (2014) is reviewed. Finally, selected numerical examples and their results from both modelling approaches are shown.