A scheme in Fig. l is used to simulate the earthquake sequence in North China in the last 700 years. With a postulated fault system, a mechanical parameter distribution and suitable boundary tractions, the stress field is computed by finite element method in plane strain. The fault zones are treated as elastic-perfectly plastic medium which flows according to Coulomb’s criterion and its associated flow rule. An earthquake is simulated by lowering the static coefficient of friction to a kinetic one at the epicenter to cause a slip there and arrived at a new equilibrium position. The differences in displacement, stress, etc. before and after the lowering are considered to be due to the earthquake and compared to observations. After adjustments of parameters and boundary tractions, the results are satisfactorily compared with the major earthquake sequence of 14 events M_s≥7 and regions of seismic risk in the last 700 years. Possible localities of future earthquake risk after the Tangshan earthquake and ef A scheme in Fig. L is used to simulate the earthquake sequence in North China in the last 700 years. With a postulated fault system, a mechanical parameter distribution and suitable boundary tractions, the stress field is computed by finite element method in plane strain. The fault zones are treated as elastic-perfectly plastic medium which flows according to Coulomb’s criterion and its associated flow rule. An earthquake is simulated by lowering the static coefficient of friction to a kinetic one at the epicenter to cause a slip there and arrived at a the differences in displacement, stress, etc. before and after the lowering are considered to be due to the earthquake and comparable to observations. After adjustments of parameters and boundary tractions, the results are satisfactorily compared with the major earthquake sequence of 14 events M_s≥7 and regions of seismic risk in the last 700 years. Possible localities of future earthquake risk after the Tangsha n earthquake and ef