Abstract: Understanding the viscoelastic structure of subduction zones is essential for assessing seismic hazards and understanding subduction-zone dynamics. However, the influence of lateral variations in elastic upper-plate thickness (H_c) remains poorly constrained and is often overlooked. In this study, we use two-dimensional forward viscoelastic earthquake-cycle models to fit both horizontal and vertical Global Navigation Satellite System (GNSS) observations. We identify a clear trade-off between locking depth (D) and H_c in both components. To resolve this ambiguity, we incorporate constraints from thermal models and tremor distributions along the Cascadia Subduction Zone. As a novel result extending beyond previous kinematic models, our results reveal a systematic northward increase in H_c from ~20 km to ~30 km. This trend correlates with increasing oceanic plate age and likely reflects variations in the subaccretion and wedge-cooling processes along the trench-parallel direction. In contrast, D remains relatively uniform at ~10 km, consistent with previous findings. These results demonstrate the robustness of our approach for simultaneously constraining H_c and D, and they suggest it may be applied to other subduction zones. Lateral variations in H_c significantly affect crust deformation and should not be ignored in earthquake-cycle models. Accounting for these heterogeneities improves estimates of H_c and D and enhances our understanding of megathrust locking, seismic hazard potential, and the physical conditions controlling episodic tremor and slip events.