## 2 Materials and Methods
### Keyboard Model of Seismic Cycles of Great Subduction-Associated Earthquakes
(**a**) Model input parameters; (**b**) block thickness; (**c**) 1-D elastic model; (**d**) time-dependent viscosity of the contact layer; (**e**) stresses acting on lateral sides of the block; (**f**) stress acting at the block’s bottom; (**g**) stress equilibrium equation; (**h**) equation governing displacement evolution; (**i**) boundary conditions applied at the frontal (\(x=0\)) and rear (\(x=1\)) edges of the block; (**j**) calculating incremental displacements; (**k**) converting displacements into elastic energy accumulated by the block; (**l**) condition to identify the instant of the earthquake; (**m**) stress dumping conditions; (**n**) output dataset. Mechanically, the mean compressional stress in \(x\)-direction is proportional to x-component of deformation (Figure 2c), with \(E\) being the elastic modulus of the block. The block is also experiencing tangential stresses \(\tau_{i-1}\) and \(\tau_{i+1}\) arising due to the friction of the blocks' lateral boundaries with neighboring blocks (Figure 2e) as well as the stress \(\tau_{0}^{(i)}\) acting on the bottom surface due to the viscous contact layer, Figure 2f. The latter depends on the strain rate in the contact layer with thickness \(h\), written as \(\dot{\gamma}=\left(\frac{\partial w_{i}}{\partial t}-V_{i}\right)\frac{1}{h}\) (\(V_{i}\) is subduction rate beneath the block), and effective viscosity of the contact layer \(\eta\left(\dot{\gamma},\,x,t\right)\) (Figure 2d). The viscosity \(\eta\left(\dot{\gamma},x,t\right)\) is described by a monotonic function of the strain rate, reflecting the transition from the stationary viscosity \(\eta_{f}\) describing the rheology of the contact layer during the interseismic phase of the cycle to a lower postseismic phase viscosity \(\eta_{a}\), associated with the contact layer destruction caused by the recently-occurred earthquake. Here \(\Gamma\) denotes the viscosity recovery parameter (\(\Gamma=10^{-4}\)\(\frac{1}{y^{\prime}}\)), \(\kappa=0.2\)-0.5 yr is characteristic time required for the contact layer defects to be healed, and \(T_{\text{ps}}\) is the time passed since the last earthquake in a given block. The lateral side stresses (Figure 2e) depend on the thickness of the viscous material filling the interblock faults \(h_{g}\), as well as \(\mu_{i-1}\) and \(\mu_{i}\), viscosities of that material between neighboring blocks number. The stress-strain state is assumed to satisfy the stress equilibrium equation (Figure 2g), where \(d_{i}\) is the block size along strike direction, parallel to the trench line. Table 2 shows typical values of the above parameters following [20]. Substitution of expressions for stresses into equilibrium equation yields the series of equations in terms of displacement \(w\) (Figure 2h). Additionally, the boundary conditions are applied (Figure 2i), which have the meaning that the outer edges of the frontal blocks are free, while the rear edges undergo elastic interaction with the fixed margin controlled by the rigidity coefficient \(k_{i}^{\prime}\). The above mathematical formalism is characterized by significant nonlinearity. To solve the problem, we apply an explicit finite-difference method yielding displacement increments \(\Delta W_{i,j}^{k}\) at nodal points \(x_{i},y_{j}\) within each block updated at time instants \(t_{k}\) (Figure 2j). Evolution of the block system is controlled by the conditions for dumping the accumulated elastic energy in the form of an abrupt decrease in compressional stresses at time instants \(t_{eq}^{i}\), corresponding to earthquakes, occurring when the block reaches the threshold energy \(\epsilon_{cr}\) (Figure 2k-n). Note that this condition is equivalent to an instantaneous change in displacement. Specific time instants of the earthquakes, \(t_{eq}^{i}\), are determined from the equality of accumulated energy to the threshold value (crit