The origin of 3D seismic heterogeneity in Precambrian lithosphere has been enigmatic, because temperature variations in old stable shields are expected to be small and seismic sensitivity to major‐element compositional variations is limited. Previous studies indicate that metasomatic alteration may significantly affect average 1‐D structure below shields. Here, we perform a grid search for 3‐D thermo‐chemical structure, including variations in alteration, to model published Rayleigh‐wave phase velocities between 20 and 160 s for the eastern part of the Archean Superior and Canadian Proterozoic Grenville Provinces. We find that, consistent with constraints from surface heatflow and xenoliths, the lithosphere is coolest (Moho heatflow 12‐17 mW/m2) and the thermal boundary layer thickest (>250 km) in the northeastern Superior and warmest in the southeastern Grenville (Moho heatflow 20‐25 mW/m2, thermal boundary thickness 160‐200 km). Compositionally, the phase velocities for most of the Superior within our study region require little alteration, but in a few regions, fast velocities need to overlie slower velocities. These can be modelled with an eclogite layer in the mid lithosphere, consistent with active‐seismic and xenolith evidence for remnants of subducted Archean crust. The phase velocities from the Grenville Province require significant metasomatic modification to explain the relatively low velocities of the shallow lithosphere, and the required intensity of alteration is highest in parts of the Grenville associated with arc accretion. Thus, composition of the northeastern Canadian Shield appears to reflect different stages and styles of craton assembly.