Vacancies strongly interact with twin boundaries and often change dramatically the properties of the ferroelastic material. However, the understanding of this behavior at an atomic-level is still deficient. Here we study vacancy diffusion processes across very large length- and time-scales using a combination of molecular dynamics and Monte Carlo simulations. We find that vacancies reduce their energy by residing at twin boundaries, kinks inside domain boundaries, and junctions between domain boundaries. Vacancies have the largest binding energy inside junctions and co-migrate with the motion of the junctions. For the weaker trapping inside twin boundaries, a “ghost line” may be generated because vacancies do not necessarily diffuse with moving boundaries and are left behind, leaving a trace of a previous position of the domain boundary. Needle twins act as channels for fast diffusion with almost one order of magnitude higher vacancy diffusivity than in bulk. The relative concentration of vacancies at twin boundaries (ρVa) is a function of the average vacancy concentration (CVa) with ρVa ∼ CVaα and α = 0.61, in contrast to that of immobile vacancies case with α = 0.4. The concentration of vacancies at twin boundaries is enriched ca. 5 times at low temperatures. With increasing temperature, the enrichment drops as the trapping potential at the twin boundaries decreases (thermal release). The distribution of energy-drop upon twin pattern evolution follows a power law. The exponent ε increases from ∼1.44 to 2.0 when the vacancy concentration increases. The power law exponent is the same in the athermal region, while Vogel-Fulcher behavior is found at high temperatures.