This study proposes a comprehensive methodology for estimating accumulative damage of bridge structures under multiple seismic excitations, in the framework of site-specific probabilistic hazard analysis. Specifically, a typical earthquake- prone region in China is chosen to perform probabilistic seismic hazard analysis (PSHA) to find the mean annual rate (MAR) of ground motion intensity at a specific level, based on which, a mass of ground motion observations is selected to construct random earthquake sequences with various number of shocks. Then, nonlinear time history analysis is implemented on the finite element (FE) model of a RC girder bridge at the site of interest, to investigate structural responses under different earthquake sequences, and to develop predictive model for cumulative damage computation, in which, a scalar seismic intensity measure (IM) is adopted and its performance in damage prediction is discussed by an experimental column. Furthermore, a mathematic model is established to calculate occurrence probability of earthquakes with various number of shocks, based on PSHA and homogenous Poisson random process, and a modified cumulative damage indicator is proposed, accounting for probabilistic occurrence of various earthquake scenarios. At end, the applicability of the proposed methodology to main shock and aftershock scenarios is validated, and characteristics of damage accumulation under different multiple earthquake scenarios are discussed.