Abstract: This paper describes the application of the Galerkin finite element method and the treatment of unknown heat flow boundary condition by hybrid method to compute the crustal temperature distribution in North China Bay area. It has been shown that: (1) temperature along the Moho discontinuity is not a constant as generally assumed, but varies from a lower value of about 700℃ at its downwarping portion to a higher value around 810℃ at its upwarping portion; and (2) large shallow earthquakes occurred along a belt of steep thermal gradient just above a low velocity (i. e. low resistivity) zone where the temperature is about 600℃, this may be due to the high temperature below such a zone, weakening the materials of lower crust unable to concentrate stress causing large earthquakes.