High-pressure metamorphic rocks provide evidence that in subduction zones material can return from depths of more than 100 km to the surface. The pressure-temperature paths recorded by these rocks are variable, mostly revealing cooling during decompression, while the time constraints are generally narrow and indicate that the exhumation rates can be on the order of plate velocities. As such, subduction cannot be considered as a single pass process; instead, return flow of a considerable portion of crustal and upper mantle material must be accounted for. Our numerical simulations provide insight into the self-organizing large-scale flow patterns and temperature field of subduction zones, primarily controlled by rheology, phase transformations, fluid budget, and heat transfer, which are all interrelated. They show the development of a subduction channel with forced return flow of low-viscosity material and progressive widening by hydration of the mantle wedge. The large-scale structures and the array of pressure-temperature paths obtained by these simulations favorably compare to the record of natural rocks and the structure of high-pressure metamorphic areas.