Seabed instability is one of the important reasons for offshore structure damage. Unlike most previous studies that treated the oscillatory and residual response separately, a coupled model for wave-induced response in non- homogeneous seabeds is proposed in the present study. Effects of spatial derivative terms in seabed parameters are introduced into the accumulation of pore pressure. Model validations are conducted by comparing the present simulation with the previous analytical solutions, wave flume tests, and numerical simulations. The validated model is applied to investigate the effects of grain size, non-homogeneous distribution of seabed parameters, and non-linear wave conditions on the wave-induced seabed dynamic response and liquefaction. It is found that (1) the oscillatory mechanism in pore pressure variation dominates in the coarser seabed, while the residual mechanism becomes obvious with the decreasing grain size, (2) consideration of the non-uniform permeability and Young's modulus would promote and suppress the pore pressure accumulation and liquefaction, respectively, and (3) the simulation error in pore pressure between homogeneous and non-homogeneous seabeds increases with the increase of the wave nonlinearity.