Abstract. This study is focused on multistable slip of earthquakes based on a one-degree-of-freedom spring-slider model in the presence of thermal-pressurized slip-weakening friction and viscosity by using the normalized equation of motion of the model. The major model parameters are the normalized characteristic displacement, U c , of the friction law and the normalized viscosity coefficient, η, between the slider and background plate. Analytic results at small slip suggest that there is a solution regime for η and γ (= 1/U c ) to make the slider slip steadily. Numerical simulations exhibit that the time variation in normalized velocity, V /V max (V max is the maximum velocity), obviously depends on U c and η. The effect on the amplitude is stronger due to η than due to U c . In the phase portrait of V /V max versus the normalized displacement, U/U max (U max is the maximum displacement), there are two fixed points. The one at large V /V max and large U/U max is not an attractor, while that at small V /V max and small U/U max can be an attractor for some values of η and U c . When U c <0.55, unstable slip does not exist. When U c ≥ 0.55, U c and η divide the solution domain into three regimes: stable, intermittent, and unstable (or chaotic) regimes. For a certain U c , the three regimes are controlled by a lower bound, η l , and an upper bound, η u , of η. The values of η l , η u , and η u − η l all decrease with increasing U c , thus suggesting that it is easier to yield unstable slip for larger U c than for smaller U c or for larger η than for smaller η. When U c <1, the Fourier spectra calculated from simulation velocity waveforms exhibit several peaks, thus suggesting the existence of nonlinear behavior of the system. When U c >1, the related Fourier spectra show only one peak, thus suggesting linear behavior of the system.