Plate tectonics is a key feature of the dynamics of the Earth’s mantle. By taking into account the stress-history-dependent rheology of mantle materials, we succeeded in realistically producing tectonic plates in our numerical model of mantle convection in a three-dimensional rectangular box. The calculated lithosphere is separated into several pieces (tectonic plates) that rigidly move. Deformation of the lithosphere caused by the relative motion of adjacent plates is concentrated in narrow bands (plate margins) where the viscosity is substantially reduced. The plate margins develop when the stress exceeds a threshold and the lithosphere is ruptured. Once formed, the plate margins persist, even after the stress is reduced below the threshold, allowing the plates to stably move over geologic time. The vertical component of vorticity takes a large value in the narrow plate margins. Secondary convection occurs beneath old tectonic plates as two-dimensional rolls with their axes aligned to the direction of plate motion. The surface heat ﬂow decreases with increasing distance from divergent plate margins (ridges) in their vicinity in the way the cooling half-space model predicts, but it tends towards a constant value away from ridges as observed for the Earth because of the heat transport by the secondary convection.