We present a new approach for the lithosphere/asthenosphere interaction in subduction zones. The lithosphere is modeled as a Maxwell viscoelastic body sinking in the viscous asthenosphere. Both domains are discretized by the Finite Element Method and we use a staggered coupling method. The interaction is provided by a non-matching interface method called the Fictitious Domain Method. We describe a simplified formulation of this numerical technique and present 2-D examples and benchmarks. We aim at studying the effect of mantle viscosity on the cyclicity of slab folding at the 660 km depth transition zone. Such cyclicity has previously been shown to occur depending on the kinematics of both the overriding and subducting plates, in analog and numerical models that approximate the 660 km depth transition zone as an impenetrable barrier. Here we applied far-field plate velocities corresponding to those of the South-American and Nazca plates at present. Our models show that the viscosity of the asthenosphere impacts on folding cyclicity and consequently on the slab's dip as well as the stress regime of the overriding plate. Values of the mantle viscosity between 3 and 5 × 1020 Pa.s are found to produce cycles similar to those reported for the Andes, which are of the order of 30-40 Myrs (based on magmatism and sedimentological records). More over, we discuss the episodic development of horizontal subduction induced by cyclic folding and, hence, propose a new explanation for episodes of flat subduction under the South-American plate.