## 1 Introduction  
Large-magnitude subduction-associated earthquakes with M > 8 cause the release of enormous elastic stresses accumulated over hundreds of years or even a millennium. Forecast of such shocks, leading to substantial social, economic and ecological consequences, is one of the most crucial problems in geophysics [1, 2]. To date, the most developed forecast-related approaches include deformational monitoring [3, 4, 5], identification of geophysical precursors, as well as statistical analysis of the seismicity pattern [6, 7, 8]. Machine learning (ML) techniques and big data analysis have also become a trend in recent years, and have potential to increase forecast reliability [9]. However, given the rare occurrence of catastrophic M > 8 events (around once-in-a-century), with only few such megaearthquakes recorded over a history of seismological observation, this approach is faced with difficulties due to a small size of the training sample. Thus, some physical (mechanical) model is strongly needed for appropriate description of geodynamical setting during the earthquake preparation phase, and such model might be employed in combination with ML methods to better constrain the results. Geomechanical models used in seismology and geodynamics are mainly based on continuum mechanics, and widely employ the concept of the continuous structure being ruptured at the instant of an earthquake which results in the formation of a plane fault [10,11]. Models of seismic cycles associated with the strongest subduction earthquakes are mostly based on the concepts of the widely used \"asperity\" model [12], which explains the difference in characteristic sizes of sources and recurrence periods of catastrophic earthquakes in various subduction zones. Later, a set of mechanical models describing the generation of strongest earthquakes with the source length up to 1000 km was proposed assuming the idea of synchronous failure of several adjacent asperities: from a first simple model [3] to a very complex model by [4,5]. An explanation of seismic-aseismic slip patterns observed on megathrust faults proposed in [13] is also worth mentioning: the results of numerical simulations based on the quasi-dynamic asperity fault model with the slip-dependent friction and stress-dependent healing quite accurately reproduces the main features of real megathrust fault behavior. However, this approach is still lacking accuracy when modeling the processes of the earthquake preparation, stress release, and post-quake stress relaxation. Analysis of the large-magnitude tsunamia earthquakes and seismotectonic features of various subduction zones suggests that these regions have essentially a fault-block structure, with certain patterns in their kinematics [14,15]. The studies based on GPS-measured displacements of the Earth's surface in the source areas of large recent earthquakes during different stages of the seismic cycle, have also confirmed the fault-block structure in the frontal part of the continental margin [16, 17, 18, 19]. Thus, geodetic observations combined with geological and seismological studies, clearly indicate that the continental lithosphere in subduction regions has complex structure, revealing a number of separate blocks which move independently from each other. Explanation of the observed displacements and velocities, which sometimes have opposite direction in neighboring blocks, is hardly possible in terms of the well-known asperity model [12]. At the same time, interpretation of this motion pattern might be possible assuming the so-called \"keyboard\" model proposed by [PERSON] et al. [20]. This concept is essentially a geomechanical model describing seismotectonic evolution of subduction regions, taking into account the fault-block structure of the continental margin. Application of the above model implies identification of large fault-separated blocks in the considered area based on geological, geophysical, seismological and geodetic data, which allows tracing their lateral boundaries.