Abstract: Slopes are highly susceptible to collapse under seismic action and are often accompanied by distinct three-dimensional features. The three-dimensional (3D) stability of unsaturated soil slopes is investigated based on the principle of limit analysis upper bound and the one-dimensional stable infiltration model. A new 3D horizontal slicing method is proposed, which effectively considers the nonlinear distribution characteristics of unsaturated soil gravity, seismic acceleration and apparent cohesion. The seismic inertia force is expressed by a modified pseudo-dynamics method (MPDM) considering soil damping and resonance. The explicit expression for the slope safety factor is derived by gravity increase method (GIM). Comparison and validation with the existing research results and a series of parameter studies were carried out. The results show that when the slope is subjected to seismic wave action close to the intrinsic frequency of the soil body, resonance phenomenon occurs and the slope safety coefficient decreases rapidly; under the action of the same seismic frequency, the safety coefficient decreases with the increase of seismic acceleration coefficient. When B/H<3, the 3D effect is obvious, and the 3D effect should be taken into account in the design of the slope; the existence of suction contributes to maintaining slope stability.