Inside subducting slabs the interaction of different metastable phase transformations (α-olivine → β-spinel, β-spinel → γ-spinel, γ-spinel → perovskite + magnesiowüstite) associated with the latent heat release and absorption for these transitions can result in a complex thermal slab structure. Thermokinetic coupling processes cause thermal anomalies and buoyancy contrasts between slab and mantle that strongly influence the stress field within subducting plates. We have described the thermal field within the downgoing slab in a Lagrangian framework by a self-consistent, high-resolution, two-dimensional thermokinetic coupling model, which also incorporates the effects of nucleation site saturation. Using an adaptive method, the threshold for the nucleation site saturation is found to be less than 0.1% transformation degree in the cold slab interior. This can significantly shift the metastable wedge of α-olivine deeper, as suggested from previous thermokinetic slab models. However, this effect is compensated by the latent heat release. During the α-olivine → β-spinel and β-spinel → γ-spinel phase transitions, latent heat release produce very sharp phase boundaries in the slab. For fast slabs the phase boundaries reveal significant metastable perturbations from the equilibrium state in the cold interior. However, the strong thermal interaction between both ongoing phase transitions results in sharp overlapping phase boundaries with dramatic consequences for the slab stress field. We have investigated the detailed physical properties and the state of stress in the deeper portion of a subducting plate by using up-to-date temperature field and physical properties determined from recent high P-T experiments. We find that the slab is denser up to 100 kg m −3 than the surrounding mantle. Remarkable denser portions up to 260 kg m −3 are located just below the uplifted α → β and β → γ equilibrium phase transitions, while lighter portions up to −210 kg m −3 can be found just above the depressed γ-spinel to perovskite + magnesiowüstite transition and for the metastable wedge of α-olivine. These density differences effectively act as a body force and produce a significant stress field, which we calculated with a two-dimensional finite-element code. The results show that the buoyancy-induced forces produce maximum shear stress up to 23 MPa along the metastable wedge and deeper portion just above the depressed last phase transition involving γ-spinel. For the latter, the dominant state of stress is down-dip compression. The calculated P-T dependent state of stress is very similar to the depth distribution of deep-focus