Twin-twin interactions significantly influence the mechanical properties of magnesium and its alloys. A thorough understanding of the underlying mechanisms governing these interactions is essential for designing Mg alloys with enhanced strength and toughness. In this work, a crystal plasticity-twin coupled phase field (CP-TPF) model incorporating multiple extension twin variants and considering the role of dislocation slipping is proposed to investigate the interactions between/among the same twin variants and those between co-zone twin variants in Mg single crystals. The model incorporates an additional energy term to represent the interaction among different twin variants and couples the CP and TPF models through order parameters and stress tensor. The simulated results show that the interaction between the same twin variants can either promote or inhibit the twin propagation, and multiple twins tend to generate concurrently in Mg single crystal to minimize the free energy associated with the accumulation of elastic strain. During co-zone twin-twin interaction, localized thickening of the recipient twin occurs due to the concentrated stresses induced by the intrusion twin, and the mutual extrusion of the two twins leads to blunting of the intrusion twin tip. Both the coalescence of the same twin variants and the formation of twin-twin boundaries between the co-zone twin variants contribute to the effective mechanism of twinning-induced hardening. Moreover, local dislocation accommodation plays a crucial role in twin-twin interactions. It relaxes the stress concentration near the twin tips and twin-twin boundaries and significantly contributes to the uneven migration of the twin boundary.