Abstract: Variations in the flow field of the surrounding of liquefiable seabed have a significant impact on the dynamic stability of subsea pipelines. The motion of seawater (waves) and the pore fluid within the seabed are described using the Reynolds-Averaged Navier-Stokes (RANS) equations and Biot's dynamic consolidation equation, respectively. Based on the extended Masing's rule, a visco-elastoplastic constitutive model capable of simulating cyclic softening and large deformation due to liquefaction of the seabed is constructed. By achieving real-time transfer of hydrodynamic pressure at the seawater-seabed interface, an integrated numerical model for wave-seabed-pipeline interaction is established. The model is validated against Sumer et al.'s wave flume experiments. Based on a typical borehole in the Jinqimen area, numerical analysis is conducted on the dynamic response of pipelines under long-term wave loading. The results show that: (1) The presence of pipelines alters the spatial distribution of liquefied zones in the seabed. Liquefaction initiates at the seabed surface and beneath the pipeline, with subsequent liquefied areas evolving upwards along the outer wall of the pipeline; (2) The rate of reduction in deviatoric stress is more pronounced in the near-field compared to the far-field, leading to significant strain concentration around the pipeline. Additionally, the cyclic shear stress ratio (CSSR) near the pipeline is higher than in the far-field, resulting in faster liquefaction rate in the near-field; (3) Intense seabed-pipeline interactions increase the depth and extent of liquefaction in the near-field seabed, increasing buoyancy forces on the pipeline while reducing frictional resistance. Plastic flow of soil within the affected zone causes soil accumulation and compression under the pipeline, collectively resulting in pipeline uplift and lateral displacement.