Short- and intermediate-wavelength gravity and geoid anomalies are used to provide constraints on the mechanical structure of subduction zones and on the forces involved. This study is based on 2D Cartesian dynamically self-consistent models with Newtonian or power law rheologies. We show that both strong decoupling of the two convergent plates (shear stresses of the order of 107 Pa) and weakened bending lithosphere are necessary to reproduce the observed geoid and gravity data. Good fits are found for relatively low failure stresses (≈30–50 MPa). For all models providing reasonable predictions of gravity, only a small fraction of the downgoing slab weight is transmitted to the surface plates. About 10% of the energy is dissipated in the contact zone between the two plates, 10% to 20% in the bending region, and more than 70% in the sublithospheric mantle. The basal tractions (on the order of ≈1−4 × 1012 N/m) induce a net motion of the plates, with the subducting lithosphere moving faster than predicted by the no-net-motion principle. A marked positive geoid anomaly is predicted above subduction zones at intermediate wavelength (λ = 2000–4000 km) in the case of pure whole mantle convection. Such large geoid highs are not observed. Introducing “partial layering” (i.e., mass anomalies hampering mantle flow through the transition zone) is necessary to reconcile model predictions and observations for these wavelengths.