[1] We use an elastoviscoplastic rheological model to study the long-term relaxation of impact crater topography on Ganymede and Callisto. We employ a complete set of rheological parameters for ductile creep in ice I h and implement a new, shallower initial shape model. Under high heat flow conditions, stresses are relieved in deeper, less viscous material, and the remainder concentrates as flexural (bending) stresses in a viscously stiffer near-surface layer. Essentially, a mechanical lithosphere forms that focuses stresses but with uplifts and stress levels that are always greater than those obtained for a purely viscous half space. We discuss how classic elastic flexure and purely viscous relaxation are special cases of this more generalized rheological model. Stress concentration in the flexing (i.e., bending) layer can exceed the plastic yield strength of ice for high enough heat flows, weakening the layer somewhat and increasing the total amount of relaxation. These plastic strains tend to be limited to small magnitudes (<0.5%) in a shallow region on the crater floor; such features may be observable in available imagery. The grain size sensitivity of the dominant creep mechanism in ice represents a new input parameter in crater relaxation studies on icy satellites (along with the elastic and plastic moduli, although to a lesser extent). Nevertheless, we find that for reasonable input parameters our simulations reproduce the continuum of relaxation states seen on Ganymede and Callisto. A combination of higher surface temperatures ($120-130 K) and larger grain sizes (>1 mm) permits both retention and relaxation of topography at the lower and higher ends, respectively, of the anticipated heat flow range. Consideration of relaxed craters in Marius Regio and of the palimpsests and penepalimpsests suggests ancient heat flows in excess of 10 mW m À2 , consistent with other estimates.