Abstract Interseismic coupling maps are important tools for estimating future earthquake hazard at subduction zones. Often interseismic coupling is estimated from surface velocities using a homogeneous elastic Earth model. However, inversions that incorporate viscoelastic mantle flow or elastic heterogeneity show that each can change the spatial distribution and magnitude of estimated coupling. We present here a Bayesian inversion of surface deformation for the boundaries of coupling zones at the Cascadia subduction zone that incorporates both elastic heterogeneity and a viscoelastic earthquake cycle model. We find that elastic heterogeneity deepens inferred coupling and worsens fit to surface deformation compared to homogeneous elastic models, but viscoelasticity partly offsets these effects and yields similar coupling and a better fit to the vertical data than homogeneous elastic models. We also find that, of all geometries and viscosities tested, downgoing plate thickness and oceanic mantle viscosity have the greatest impact on modeled displacements and coupling.