## 1 Introduction  
In China, high-speed railways are increasingly crossing high-intensity earthquake-prone areas given their wide distribution, thus influencing the running safety of trains [1]. The running safety problems of trains under earthquake excitation has become an ongoing concern. Under earthquake excitation, extreme contact behaviors such as wheel flange contact and wheel jumping occur, directly impacting the running safety of trains [2]. Therefore, it is crucial to accurately simulate the wheel-rail extreme contact state. Wheel-rail contact modeling under earthquake excitation needs to consider wheel-rail spatial rolling contact characteristics, and the typical wheel-rail contact model is associated with three problems: (i) wheel-rail contact geometry calculation, (ii) normal contact force calculation, and (iii) tangential contact force calculation [3]. In previous studies, the contact trace method [4], nonlinear Hertz theory, and Shen-Hedrick-Elkins theory [5] have been used to solve the above problems and establish a corresponding model [6; 7; 8]. Since this model can balance the requirements of wheel-rail contact modeling and computational efficiency [9; 10], it is widely applied in vehicle dynamic simulations and running safety analyses under earthquakes [11; 12; 13]. However, in order to ensure the accuracy of the wheel-rail extreme contact state simulation, the algorithm used in wheel-rail contact modelling under earthquake excitation in the following aspects needs to be further clarified. First, the contact points are generally detected by the maximum wheel-rail vertical penetration in most cases; thus, only a single contact point can be determined [14]. Seismic excitation causes a large relative lateral displacement between the wheel and rail [15], leading to simultaneous elastic penetrations in different regions [16; 17]. Hence, the consideration of the number of wheel-rail contact points under earthquake excitation needs to be further clarified. The wheel-rail contact geometry relationship, as the basis for contact force calculation, mainly includes determining the contact point and the normal compression amount, and the wheel-rail searching method plays a decisive role. The vertical searching method [9], in which the contact point and the normal compression amount are determined according to the maximum wheel-rail vertical penetration, is widely adopted for the convenience of calculation. As the normal compression amount is a key parameter in calculating the contact normal force [17], it is essential to ensuring its calculation accuracy. The approximate calculation method based on vertical penetration is applicable to ordinary conditions [18], whereas it is not convincing under extreme operation conditions such as during earthquakes. The normal searching method [19; 20; 21], in which the contact point position and the normal compression amount are determined according to the maximum wheel-rail normal penetration, can ensure the reliability of the solution to the normal compression amount under all conditions. Hence, it is crucial to clarify the necessity of introducing the normal searching method for seismic conditions by comparing the two methods. The normal contact stiffness, which is another key parameter in the calculation of the wheel-rail normal contact force, is calculated on the basis of the geometric parameters at wheel-rail contact points when using the Hertzian contact theory [22]. However, an empirical formula based on wheel-rail constants [23] is widely used in the wheel-rail contact force calculation, which is only applicable in cases where the contact points are located in a narrow region of the wheel tread and rail head [24]. Hence, the normal contact stiffness algorithm also needs to be further clarified for seismic conditions.