Empirical Analysis of Global-Scale Natural Data and Analogue Seismotectonic Modelling Data to Unravel the Seismic Behaviour of the Subduction Megathrust

Subduction megathrusts host the Earth’s greatest earthquakes as the 1960 Valdivia (Mw 9.5, Chile), the largest earthquake instrumentally recorded, and the recent 2004 Sumatra-Andaman (Mw 9.2, Indonesia), 2010 Maule (Mw 8.8, Chile), and 2011 Tohoku-Oki (Mw 9.1, Japan) earthquakes triggering devastating tsunamis and representing a major hazard to society. Unravelling the spatio-temporal pattern of these events is thus a key for seismic hazard assessment of subduction zones. This paper reviews the current state of knowledge of two research areas–empirical analysis of global-scale natural data and experimental data from an analogue seismotectonic modelling—devoted to study cause-effect relationships between subduction zone parameters and the megathrust seismogenic behavior. The combination of the two approaches overcomes the observational bias and inherent sampling limitations of geological processes (i.e., shortness of instrumental and historical data, decreasing completeness and resolution with time into the past) and allows drawing appropriately from multiple disciplines with the aim of highlighting the geodynamic conditions that may favor the occurrence of giant megathrust earthquakes.


INTRODUCTION
Most of the global seismicity and the largest earthquakes occur at subduction zones where oceanic lithosphere is recycled as it sinks into the underlying mantle. The downgoing movement of the subducting plate and friction acting at the shallow part of the interface between the subducting and overriding plate (aka megathrust) cause shear stresses to accumulate slowly over centuries. When the frictional shear strength of the interface is overcome, these stresses are released episodically by a variety of seismic and aseismic slip modes, including "mega-earthquakes": seismic megathrust rupture causing earthquakes with magnitudes (M W ) ≥ 8.5. Stress accumulation and release occurs in cycles known as seismic cycles. In the last 15 years, a cluster of mega-earthquakes occurred at an increased rate with respect to the last century (e.g., Ammon et al., 2010) often revealing unexpected characteristics (e.g., Stein and Okal 2011;Yue and Lay 2013) and resulting in catastrophic direct but also indirect (e.g., tsunamis) effects as in Sumatra (2004) (e.g., Subaraya et al., 2006), and Tohoku-Oki 2011 (e.g., Romano et al., 2012). The attention for these events has thus recently increased, raising interest on factors controlling their spatio-temporal pattern, including their potential maximum magnitude and return periods which are fundamental parameters in earthquake and tsunami hazard assessment. This is a challenging goal not only because each subduction zone has unique characteristics (e.g., Kopp, 2013) but also because the seismogenic portion of the megathrust is located offshore, below the sea-level , and direct observation is restrained to shallow drilling (e.g., Japan Trench Fast Drilling Project; Chester et al., 2012) and the study of exhumed prototypes (e.g., Remitti et al., 2011). Moreover, the instrumental seismic record is limited to the last century, a very short time interval compared to the recurrence time of great earthquakes, which ranges from centuries to millennia (e.g., McCaffrey, 2008). In some cases, seismic records of megathrust ruptures can be extended back in time with paleoseismological investigations (e.g., Philibosian and Meltzner 2020), however such data typically lack the spatio-temporal completeness and resolution needed for rigorous statistical analysis.
One of the empirical research strategies that allow overcoming limited records consists in trading time with space at the global scale: given that different subduction zones are at individual stages of their seismic cycle, we may be able to infer the general seismogenic behavior of subduction zones by integrating information from worldwide seismicity (e.g., Schafer and Wenzel, 2019). The second approach, somehow complementary to the first one, is comparing the seismic character of individual megathrusts with the widest range of geological, geophysical and geodynamical parameters, with the aim of statistically identifying the combination of factors controlling subduction seismicity (e.g., Ruff and Kanamori, 1980;Peterson and Seno, 1984;Heuret et al., 2011;Schellart and Rawlinson, 2013;Brizzi et al., 2018).
Results from such a global-scale natural data analysis approach can be then tested against experimental and modelling data in order to validate inferred causal relationships by means of ad hoc parametric studies. Analytical (e.g., Ruff, 1992), numerical (e.g., van Dinther et al., 2013, physical (e.g., Den Hartog and Spiers 2013) and analog models (e.g., Rosenau et al., 2017) can reproduce a long series of seismic cycles and, by tuning ad hoc boundary conditions, they allow testing the role played on seismicity by the different parameters deduced from the analysis of natural data (e.g., Corbi et al., 2017b). Despite models are only simplified representations of the complexity of the natural prototype due to the unavoidable computational and experimental limitations, their results are inspiring for the interpretation of the limited and sparse natural data, thus enhancing our current understanding of subduction seismicity.
We focus our brief overview on the state of the art of the analysis from global-scale natural data and experimental data from analogue seismotectonic modelling, highlighting the advantages of their use and the great potential for future developments in the study of the earthquake process that can be drawn from the close interaction and reciprocal feedback between the two methods. Very recently, also numerical methods have been developed to setup simulations which capture both the long-term and the short-term evolution of the subduction process. While numerical simulations proved particularly useful for testing hypothesis and mapping a wide parameter space, it is the emergent behavior of analogue models which makes them appealing for an empirical analysis as addressed in this review.

OBSERVATIONS AND INFERENCES FROM NATURAL DATA
There is a noteworthy variability in the maximum moment magnitude of earthquakes originated at the subduction megathrust during the instrumental era (e.g., Dixon and Moore, 2007 Figure 1A): some subduction zones-e.g., Alaska, Chile, Sumatra, NE Japan-experienced mega-earthquakes, while others-e.g., Tonga Mariana-did not ( Figure 1A). What is unclear and still debated is if this dichotomy results only from the short observational time span and therefore all the worldwide subduction zones could host mega-earthquakes given a long enough observation period (e.g., McCaffrey, 2008) or if specific local conditions are responsible for the diversity of the interplate seismicity (e.g., Ruff and Kanamori, 1980). It is thus tempting to identify which geological/geophysical features, individually or in combination, possibly have a cause-effect relationship with megathrust seismicity (e.g., Pacheco et al., 1993;Heuret et al., 2011;Schellart and Rawlinson, 2013) and to quantify the likelihood for great earthquakes along worldwide subduction megathrusts (e.g., Marzocchi et al., 2016).
The "comparative subductology" was the intuition to disentangle the overall behavior of subduction by comparing a wide range of observables (i.e., subducting plate age and dip, convergent velocity, backarc tectonic style, characteristics of arc volcanism and seismicity, accretionary/erosional character of the margin) with seismicity characteristics of convergent margins (Uyeda, 1982) with the aim of deducing possible causal relationships ( Figure 1B). This space-time trade-off approach allowed to overcome limitations due to uniqueness of local conditions (e.g., Suárez and Albini, 2009;Becker and Meier, 2010;Hayes and Furlong, 2010) and to the undersampling of the modern seismicity record.
Early on, two end-members subduction types -Mariana and Chileanhave been proposed, showing also significant difference in the degree of mechanical coupling (i.e., the ability to accumulate stresses) along the plate interface and in the consequent capability to host mega-earthquakes (Uyeda and Kanamori, 1979). The Chilean-type boundaries seemed more coupled and prone to trigger mega-earthquakes than the rather aseismic Marianas, with the different behaviors likely tuned by both different densities of the subducting plate and motions of the upper plate (Uyeda and Kanamori, 1979).
Based on this idea, Ruff and Kanamori (1980) looked at the maximum M w recorded in the instrumental era and they found a striking correlation between the convergence rate and plate age. It was thus proposed that mega-earthquakes preferentially occur where a young lithosphere is converging fast with the overriding plate, allowing for a high mechanical plate coupling along the megathrust. This relationship has long been supported by studies involving larger, progressively updated and refined datasets (e.g., Peterson and Seno, 1984;Jarrard, 1986;Pacheco et al., 1993;Scholz and Campos, 1995). However, variations of coupling found within individual subduction zones (Scholz and Campos, 2012), the identification of different moment release rates in subduction zones sharing similar lithospheric ages and speeds (Peterson and Seno, 1984) and, last but not least, updated data (e.g., Heuret et al., 2011;Stein and Okal, 2007;Stein and Okal, 2011) initiated skepticism in the acknowledgement of the Kanamori (1980) model (e.g., McCaffrey, 1994) that was finally challenged by recent mega-earthquakes occurring in locations assumed to be weakly coupled (e.g., Wang, 2013). Among these, the Tohoku-Oki earthquake stands out not only for exceeding the regional maximum magnitude earthquake for that region according to Ruff and Kanamori (1980) model, but also as case of extreme violation of commonly accepted scaling laws (e.g. Okal, 2015).
A new phase of the comparative subductology involved not only the use of large dataset (e.g., Heuret et al., 2012;Schellart and Heuret et al., 2011). M w ≥ 8.0 subduction interface events are plotted by stars of different colours, representing the interval of occurrence of the events. The largest events are not homogeneously distributed in space and time along the worldwide subduction zones; (B) Correlation matrix used to show the dependence between multiple variables at a glance. The table contains the modulus of the Pearson's product moment correlation coefficient R of linear bivariate regressions between seismogenic zone parameters (t_tot: total seismic of subduction zones; t_inter: seismic rate due to interplate earthquakes; Mmrr_tot: total moment release rate of subduction zones; Mmrr_inter: moment release rate due to interplate events; Mmax_inter: maximum magnitude of interplate events; Cn: seismic coupling coefficient; zmax_inter: downdip limit of the subduction thrust fault; zmin_inter: updip limit of the subduction thrust fault) and other independent subduction parameters (Tsed: thickness of sediments at the trench; UpPlate: nature of the upper plate; thermal: thermal parameter; age: age of the subducting lithosphere at the trench; Vsn: normal component of the subduction velocity; Vsupn: normal component of the overriding plate velocity; Vts: normal component of the trench velocity; Vsubn: normal component of the subducting plate velocity; Oc: convergence obliquity; Os: subduction obliquity; R: radius of bending of the slab) as attempt to unravel the overall behavior of subduction by comparing a wide range of observables (see Heuret et al., 2011 for additional details).

INSIGHTS FROM ANALOG MODELLING OF SUBDUCTION MEGATHRUST SEISMIC CYCLES
To overcome limitations of natural records and test cause-effect relationships suggested by global-scale data analysis we promote an analog modelling approach. Analog models are down-scaled and simplified versions of natural subduction zones used in experiments that allow reproducing the essence of a variety of geological phenomena, in convenient spatial and temporal scales (e.g., Hubbert, 1937;Funiciello and Corbi, 2017). Despite analog models might appear oversimplified with respect to high performance computational modelling strategies anticipated in the introduction (for example, no efficient control on temperature and water content has been implemented so far), they have the advantage of evolving spontaneously on the basis of the physics of the process, which can be easy to tune and to interpret. Moreover, analog models are intrinsically threedimensional and occur in a space-time continuum, both keys in allowing across-scale processes to emerge, as e.g. the growth of mega-earthquakes (e.g., Brizzi et al., 2018). The recent advances in characterization of rheological (e.g., Di Giuseppe et al., 2014;Rudolf et al., 2016) and frictional properties (Klinkmüller et al., 2016) of analog materials and the improvement in monitoring techniques  supported the development of socalled "seismotectonic" analog models (Rosenau et al., 2017). Those models reproduce the elastic loading and sudden release in a cyclic fashion similarly to the seismic cycle of subduction megathrusts ( Figure 2). All seismotectonic analog models developed so far share the common characteristics listed below: -Kinematic boundary condition: deformation of the overriding plate is imposed by underthrusting of the subducting plate at a constant rate. -Wedge geometry: elastic stresses are stored into the analog overriding plate which might be elastic or elastoplastic while the subducting plate is assumed being rigid; the initial wedge surface and the analog megathrust are flat (i.e., with no downdip-and along trench curvature). -Rate-state friction: the analog megathrust includes both (aseismic) velocity-strengthening and (seismogenic) velocityweakening frictional regions. -Spontaneous earthquake nucleation: slip episodes occur spontaneously when elastic stress overcomes the frictional strength of the megathrust locally (i.e., there is no external triggering). -Laboratory geodetic monitoring: models are monitored via cross-correlation of images of the top or side surfaces or the interior of the model (e.g., with the Particle Image Velocimetry, PIV method; Sveen 2004). Applied to the top surface this technique is equivalent to having a dense network of continuous GPS stations (hence the name "laboratory geodesy models") which allow inverting for relevant earthquake source parameters (e.g., Rosenau et al., 2019;Corbi et al., 2013;Kosari et al., 2020).
The first articles about analog modelling of megathrust seismic cycling focused on reporting rheological and frictional properties of analog materials (Corbi et al., 2011) and their implementation into a scale model . The experimental behavior was composed of the archetypical pattern of alternating relatively longer periods of low, landward deformation rates with sudden, higher, trenchward velocities reversals as in the seismic cycle of subduction megathrusts Corbi et al., 2013; Figure 2 panels c-f). Follow-up studies focused on applications demonstrating the capabilities of analog models in studying a variety of earthquake phenomena, including recurrence pattern Rosenau et al., 2019) and the moment-duration proportionality . Analog models have been used for investigating extrinsic parameter control, for example how the width of the seismogenic zone (Corbi et al., 2017a), the upper plate rheology (Brizzi et al., 2016), the subduction velocity (Corbi et al., 2017a) or the subduction megathrust roughness (van Rijsingen et al., 2019) influence the seismic behavior. Analog models helped in constraining the intrinsic variability of earthquake ruptures and associated tsunamis, validating empirical tsunami forecast Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 600152 models (Rosenau et al., 2010). Analog models have been also used for identifying the geometrical and mechanical conditions leading to synchronization of earthquakes (Corbi et al., 2017b;Rosenau et al., 2019). The process of earthquake synchronization is one of the key aspects in seismogenesis and is considered responsible for the origin of megaearthquakes unzipping large parts of a subduction zone. Nowadays, analog models are also used for testing earthquake predictability (e.g., Corbi et al., 2019) and for identifying optimal geodetic configurations for monitoring convergent margins Kosari et al., 2020).

CHALLENGES AND FUTURE DIRECTIONS
Over the last decades, a wide range of studies have shed light on the seismogenic behavior of subduction megathrusts. This short review has recapped the contributions provided by the comparative analysis of global-scale natural data and analogue modelling.
Forecasting which subduction zone will host the next megaearthquake and its timing remains a major challenge. The efforts to highlight which are the conditions enhancing the "productivity" of mega-earthquakes have recognized the potential role of several subduction-related parameters, but no consensus has been reached so far. For instance, empirical results point out at different parameters such as plates curvature (Bletery et al., 2016), roughness (e.g., van Rijsingen et al., 2019) or trench-parallel length of the subduction zone (Brizzi et al., 2018) and amount of trench sediments (Brizzi et al., 2020). Convergent margins are complex tectonic settings, whose diversity results from the interplay of several parameters operating at different spatial and temporal scales and whose relative importance is thus hard to assess (e.g., Figure 1B). Moreover, the understanding of subduction megathrusts is biased by the short instrumental seismic record and by the lack of completeness and of a homogeneous resolution of the observables characterizing subduction zones.
Is it thus really feasible to unravel the seismic behavior of the subduction megathrust or is it an impossible mission?  Rosenau et al., 2017;Corbi et al., 2017a, respectively). Both setups are made of a glass-sided, wedge shaped box whose base is inclined 10°-15°in agreement with natural observables of subduction mega thrusts . Both are equipped with a basal conveyer plate moved by a motor-driven piston with variable velocities of the order of 10 4 -10 −5 m/s as analogue, analog of the downgoing subducting plate. HetTec models use a mix of granular material composed by rice, sugar and rubber pellets to simulate an elastoplastic continental lithosphere and silicone oil to simulate the mantle. LET models are made of a visco-elastic gelatin which is deformed on the velocity-weakening seismogenic zone of the subduction thrust which simulated by sandpaper at the gelatin-basal plate interface. Both models experience stick-slip dynamics [orange time series in (C)] as response to subducting plate underthrusting. Particle Image Velocimetry allows monitoring the evolution of the displacement field as in a dense, homogeneously distributed geodetic network [gray dots in (D)], spanning also above the generally offshore seismogenic zone [cyan rectangle in ( The continuous scientific and technological advancements allow detecting transient deformation signals occurring on short times-scales (fractions of the seismic cycle) along convergent margins and sampling the megathrust seismic cycle with progressively denser and higher quality data. For example, even if the underlying source process still has to be identified, both the 2010 Maule (Chile) and 2011 Tohoku-Oki (Japan) earthquakes have been preceded by potentially diagnostic geodetic transients lasting several months and spanning thousands of kilometers . Satellite measurements (e.g., InSAR studies; Jolivet et al., 2020), dense geodetic (e.g., Michel et al., 2018) and seismic (e.g., Brodsky and Lay, 2014) monitoring -both onshore and offshore-is able to record a wide range of deformation signal frequencies (i.e., from earthquakes s.s. to low frequency earthquakes and tremors) and the improvement of the quantity, quality and accessibility of databases of convergent margin observables (e.g., geodesy.unr.edu/NGLStationPages/gpsnetmap/ GPSNetMap.html; geofon.gfz-potsdam.de; www-solid.eps.s. u-tokyo.ac.jp/∼sloweq) yield a promising perspective to successfully elucidate the spatio-temporal pattern of earthquake occurrence along the worldwide subduction zones and to mitigate hazard related to mega-earthquakes and tsunamis.
Another support to better understand subduction earthquakes and to advance in the forecast of earthquakes is coming from recent improvements in data analysis methods. Machine learning is an application of artificial intelligence able to automatically discover hidden patterns learning directly from analysis of large datasets, without being explicitly programmed (e.g., Bergen et al., 2019;Ren et al., 2020). This approach already succeeded for the prediction of timing (Rouet-Leduc et al., 2018), timing and locations  and timing and duration (Hulbert et al., 2019) of slip episodes in laboratory experiments mimicking the basic physics of the seismic cycle.
The future is challenging also on the analogue modelling side, where the use of higher resolution and faster cameras and of optimized procedures for calculating surface velocities  permits the model monitoring with increasing detail. Thus, the detection threshold of analog earthquakes as well as the temporal resolution of the nucleation and rupture process are progressively improved. Laboratory geodetic methods will be accompanied by laboratory seismology using microelectromechanical sensors (MEMS) measuring model accelerations with kHz sampling rates, corresponding to the frequency range of interest for earthquakes. This will soon put us in a position to detect and localize analog earthquakes which are too small to produce a surface deformation signal measurable with optical methods. Apart from such kinematic monitoring methods, dynamic monitoring of (boundary) stresses (Ritter et al., 2018) will help in detecting events but also in understanding the relationship between stress and strain during seismic cycles. Moreover, the joint use of analogue and numerical modelling is providing successful and promising results to shed light on the behavior of the subduction megathrust , Corbi et al., 2017avan Dinther et al., 2013). The combination of these two complementary approaches allows to valorize the respective strengths: analog models are physically self-consistent, and stresses/strain evolve spontaneously in response to the applied boundary conditions, while numerical models are more adaptable and effective for parametric studies.
With respect to great subduction megathrust earthquakes, our current understanding is that they show a complex but in parts deterministic pattern. Whether (mega-) earthquakes are predictable or not can be thus answered with some hopefulness in the future. Since the seismic behavior of the subduction megathrust depends on a wide combination of parameters, interdisciplinary research needs to be strengthened integrating geological, geophysical, and geodynamical constraints and contextualizing them in a physicsbased modelling framework allowing to analyze and forecast megaearthquake triggering processes in subduction zones.
Multidisciplinary studies are always challenging, but they represent the key for a substantial advancement in the understanding of this scientifically intriguing and socially hazardous issue. In this frame, collecting new data, sharing them on open platforms, networking resources and events are fundamental steps for exchanging ideas and fostering discussion towards the comprehension of the processes controlling the megathrust seismic behavior.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.