When evaluating the seismic reliability of urban WDSs, most of traditional methods use the same approach to simulate the leakage and burst, which leads to inaccurate results of simulation and reliability evaluation. Considering pipe damages due to earthquake are random, this paper proposed a stochastic simulation method for seismic reliability assessment of WDSs. Specifically, the seismic damage scenarios are generated using Monte Carlo simulation, in which Poisson and uniform random number are used for determining pipe working status, and normal random numbers are used for pipe leakage coefficient generation. Pipe leakage and burst are simulated based on the emitter and pipe close function in EPANET. To improve hydraulic simulation accuracy of the damaged WDSs, the long pipes are segmented and "two-step iteration" method is introduced to solve the network in low pressure condition. The ratio of the statistical average values to normal values of nodal demand in seismic damage scenarios is used as reliability index. The proposed algorithm and GIRAFFE software are used to evaluate the performance of WDSs of an ancient town in Yunnan under VⅢ and IV seismic intensity, and the feasibility of the proposed algorithm is demonstrated by analyzing and comparing the evaluating results. Evaluation results indicate that, in addition to seismic intensity factor, the hydraulic condition has great influence on the seismic reliability of WDS. Under the same seismic intensity, the water supply reliability for areas located in high terrain or far from water resource is considerably lower than that of areas closed to the water source with main pipe being serviced.