Viscoelastic effects potentially play an important role during all phases of the earthquake cycle in subduction zones. However, most current models neglect such effects in the interseismic deformation pattern. Here we use finite element method (FEM) models to investigate the control of viscoelasticity on interseismic deformation and to highlight the pitfalls of interpreting the data with purely elastic models for both the forward and inverse problems. Our results confirm that elastic models are prone to overestimating the interseismic locking depth, a crucial parameter for estimating the maximum possible earthquake magnitude. The application of the viscoelastic model improves the fit to the interseismic deformation, especially in the inland area. Additionally, we construct 3‐D FEM models constrained by geophysical and GPS data and apply our methodology to the Peru‐North Chile subduction zone. Our results indicate that viscoelastic effects contribute significantly to the observed GPS data. The signals interpreted as back‐arc shortening in the elastic model can be alternatively explained by viscoelastic deformation, which, in turn, dramatically refines the interseismic locking pattern in both dip and strike directions. Our viscoelastic locking map exhibits excellent correlation with the slip distributions of previous earthquakes, especially the recent 2014 Mw 8.1 Iquique earthquake. The incorrect elastic assumptions affect the analysis of interseismic deformation with respect to slip deficit calculations. Our results thus suggest that it is necessary to thoroughly reevaluate existing locking models that are based on purely elastic models, some of which attribute viscoelastic deformation to different sources such as microplate sliver motions.

In a viscoelastic Earth, stresses slowly built up due to fault locking are relaxed concurrently during the entire interseismic period. This interseismic stress relaxation causes crustal deformation much farther away from the locked fault than can be explained using elastic models that neglect the relaxation. Here we develop a viscoelastic geodetic inversion model to address this problem at Cascadia. We invert ~500 horizontal velocity vectors based on continuous and campaign geodetic measurements over the past two decades. Ambiguities arising from long‐term rotation of upper‐plate crustal blocks are addressed by test‐correcting the geodetic velocities with two different block‐motion models. Fault back slip (i.e., slip deficit) Green's functions are derived using a Maxwell viscoelastic finite element model with realistic subduction zone structure and megathrust geometry. The preferred model features a narrow and shallow megathrust locked zone, consistent with earlier thermorheological reasoning. For an elastic model to fit the data to the same fidelity, megathrust locking has to extend to much greater depths. However, even with the viscoelastic model, the land‐based geodetic data still cannot resolve whether there is some creep (incomplete locking) in the shallowest part of the megathrust far offshore. Neither can the land data fully resolve along‐strike variations of the locking state. These ambiguities can be resolved only when adequate seafloor geodetic data are obtained.