Recent fast developments of satellite magnetic observations facilitate global Lithospheric Magnetic Field (LMF) modelling and their applications to subsurface tectonics. Here, the vertical component (Bz) of LMF at an altitude of 200km in Mainland China and surroundings is calculated from two global LMF models NGDC-720 and EMM2017. Next, Bz is used to invert the Curie Point Depth (CPD) by Equivalent Source Dipole (ESD) forward and Nonlinear Conjugate Gradient Method (NCGM) inversion scheme. Then, the surficial Heat Flux (HF) is derived by a simple one-dimensional steady heat conduction equation from the CPD distribution. At last, the continental seismicity is compared statistically to Bz, CPD and HF. Our essential conclusions are as follow: 1) Histograms and boxplots show that most (81.8%) earthquakes (EQs, Ms≥5.0) occurred in negative Bz areas, and more than a half (53.2%) number of EQs (corresponding to an energy percent of 94.6%) occurred inside areas with Bz between -5 and -3nT, in a period between 2004 and 2007, which is the same with the satellite data collection. When the time span is extended (most to 110 years), these phenomena maintain while weaken; 2) Most (88%) EQs occurred in areas with CPD between 10 and 30km, while only a few (7% and 5%) occurred in shallow (<10km) and deep (>30km) CPD areas, in a period between 2000 and 2010; 3) EQs seldom occurred inside cold areas (HF<50mW/m2), and are prone to occur in warm areas (HF>120mW/m2). EQs are also prone to occur along the boundaries of warm or cold areas. The mechanism of the correlations between EQs and Bz, CPD and HF maybe the lithospheric strength jumps caused by the temperature variations at boundaries between blocks with different CPDs.