We present a three-dimensional numerical model of plate tectonics (Figure 1; grayscale shows viscosity and arrow shows velocity). Tectonic plates steadily move on the Earth, because the lithosphere “remembers” the plate boundaries formed by its rupture that took place in the past: Once the stress in the lithosphere exceeds a threshold, the lithosphere is rifted, a sharp plate boundary develops, and the adjacent plates begin to rigidly move with respect to each other. This plate boundary and the resulting relative plate motion persist, even after the stress is reduced below the threshold around the plate boundary. We developed a model of this stress-history dependent nature of the rheology, and implemented the model to our three-dimensional numerical model of mantle convection (ACuTEMAN) in an internally heated rectangular box of 3000 km in depth and 12000 km×12000 km in horizontal dimensions. In our models, the lithosphere is rifted into several pieces, each of which rigidly moves and rotates, as observed for the Earth. The deformation of the lithosphere concentrates to narrow plate margins where the vertical component of rotation of the velocity field takes large values. The plate margins persist and allow the relative plate motion across the margin for billions of years. The heat flow takes a high value at spreading centers, and declines with the distance from the centers in their vicinity. Away from the centers, however, the heat flow tends to a steady value that is determined from the secondary convection that occurs beneath the tectonic plates as rolls with their axis aligned in the direction of plate motions.