Abstract: A new method for seismic vulnerability analysis of multi-story and multi-span underground structures based on nonlinear incremental dynamic analysis (IDA) was presented. A two-dimensional finite element model was established for a typical three-storey and three-span subway station by considering nonlinear soil-structure interaction. The equivalent linear model was adopted in the analysis to simulate the shear modulus degradation and damping characteristics of the soil under seismic excitations. The hysteretic behavior of the multistory underground structures under ground shaking was simulated by using fiber beam-column elements. A set of 21 ground motion records selected from the PEER strong earthquake record database were uniformly scaled based on the median peak ground acceleration to different intensity levels as the input ground motions at the bedrock of nonlinear IDA. It is found that the peak acceleration and peak velocity at the ground surface are efficient and appropriate intensity measures of the ground motions for the prediction of seismic response of the shallowly buried multistory underground structures and can be used to construct the seismic fragility curves of underground structures. Besides, the seismic fragility curves obtained from this numerical study were validated against the existing empirical and numerical seismic fragility functions of buried rectangular underground structures. The validated fragility curves can be used as an effective tool to quickly assess the seismic performance of the underground structures.