Rectangular underground structures in liquefiable deposits are susceptible to larger deformation during earthquake events. The intensity of earthquake motions, as one of the critical factors that would induce large deformation of underground structures in liquefiable deposits, however, was not received enough recognitions. In this paper, nonlinear dynamic analyses with a probabilistic seismic demand model (PSDM) were conducted for a two-dimensional liquefied soil-structure system, considering a variable thickness of liquefiable layer below the structure. The nonlinear soil behavior was simulated by the nonlinear constitutive model PDMY02. The constitutive parameters of saturated sand were calibrated by a slope model in the centrifuge test from LEAP-2017. Interstory drift ratio (IDR) and uplift ratio (UR) of underground structures were selected as engineering demand parameters (EDP). There were totally 18 recorded earthquake motions and 8 intensity measures (IMs) included for the PSDM. Numerical results show that the uplift behavior of underground structures is more sensitive to the variable thickness of liquefiable deposits compared to the horizontal drift of the structure. According to the performance of correlation, efficiency, practicality and proficiency, the velocity-related IMs are more appropriate than the acceleration-related IMs for the evaluation of seismic deformation behavior of underground structures. The deformation response as well as the adaptability of velocity-related IMs are further discussed and evaluated by considering more general cases including different buried depths and height-width ratios for rectangular underground structures.