What Controls Maximum Magnitudes of Giant Subduction Earthquakes?

Giant earthquakes with magnitudes above 8.5 occur only in subduction zones. Despite the developments made in observing large subduction zone earthquakes with geophysical instruments, the factors controlling the maximum size of these earthquakes are still poorly understood. Previous studies have suggested the importance of slab shape, roughness of the plate interface contact, state of the strain in the upper plate, thickness of sediments filling the trenches, and subduction rate. Here, we present 2‐D cross‐scale numerical models of seismic cycles for subduction zones with various geometries, subduction channel friction configurations, and subduction rates. We found that low‐angle subduction and thick sediments in the subduction channel are the necessary conditions for generating giant earthquakes, while the subduction rate has a negligible effect. We suggest that these key parameters determine the maximum magnitude of a subduction earthquake by controlling the seismogenic zone width and smoothness of the subduction interface. This interpretation supports previous studies that are based upon observations and scaling laws. Our modeling results also suggest that low static friction in the sediment‐filled subduction channel results in neutral or moderate compressive deformation in the overriding plate for low‐angle subduction zones hosting giant earthquakes. These modeling results agree well with observations for the largest earthquakes. Based on our models we predict maximum magnitudes of subduction earthquakes worldwide, demonstrating the fit to magnitudes of all giant earthquakes of the 20th and 21st centuries and good agreement with the predictions based on statistical analyses of observations.


Introduction and Key Observations
There are different hypotheses about the role of subduction zone parameters in producing giant earthquakes (GEQs). Based on observations available in the 1970s, Kelleher et al. (1974) proposed that the key factor controlling the magnitude of events is the width of the seismogenic zone. They suggested that slabs dipping with low angles have a larger contact area with the overriding plates, which allows for the occurrence of GEQs. Later, Ruff and Kanamori (1980) proposed that the maximum magnitude of an earthquake is controlled by two parameters: the age of the subducting plate and the plate convergence rate. They based their analysis on the sparse record of GEQs available at that time and previous works, including the study by Kelleher et al. (1974), and the analysis of the different types of subduction in the Pacific by Uyeda and Kanamori (1979). The main assumption of Ruff and Kanamori (1980) was that the magnitude of earthquakes is directly related to the strength of mechanical coupling between subducting and overriding plates. They argued that since the strength of mechanical coupling depends on the area of contact and the normal component of the stress at the fault, a combination of a relatively high rate of loading and a young age of the subducting plate (and hence low negative buoyancy) creates a wide zone with the necessary high stress on the fault plane. The view by Ruff and Kanamori (1980) became the "classical" one until the suggested correlation was violated by a number of GEQs such as the two largest earthquakes of the 21st century, the Great Sumatra/Andaman 2004 earthquake and the Tohoku 2011 earthquake (Heuret et al., 2011).
Simplified numerical modeling by Herrendörfer et al. (2015) suggested that the largest earthquakes may be related to the widest seismogenic zones. Analogue and simplified numerical modeling studies of megathrust earthquakes (Corbi et al., 2017) tested the sensitivity of magnitudes of earthquakes in asperities of different sizes to different convergence rates. That study showed that the magnitude of an earthquake is directly related to the width of the seismogenic zone and that there is no relation between subduction rate and maximum magnitude. Based on observations and analytical estimations, Bletery et al. (2016) suggested that in addition to low dipping angle, an enough small curvature of the subduction interface was a necessary condition for the GEQ's occurrence.
1.1. GEQs, Thickness of Trench Sediments, and Upper Plate Strain Ruff (1989) extended the concept of Ruff and Kanamori (1980) by suggesting that in addition to subduction velocity and slab age, the thickness of the sedimentary layer at the trench is an important factor controlling the occurrence of GEQs. He noticed that most of the giant (Mw > 8.5) earthquakes occurred in subduction zones with a thick sedimentary cover at the trench and suggested that the subduction of these sediments creates a homogeneous contact zone between the plates. Discussion of the effects of sediments in subduction zones in relation to great and GEQs has been updated with additional data (Heuret et al., 2012, and references therein), which show that earthquakes with Mw ≥ 8.5 are observed only in subduction zones where sediment thickness at the trench is larger than 0.5 km. This outcome is also supported by the observation that all subduction zones with historical earthquakes exceeding magnitude 8.5 have a rather smooth subducting sea floor (Wang & Bilek, 2014). It is inferred that thick sediments form a laterally homogeneous layer that smooths the pattern of asperities on top of the downgoing oceanic plate, thereby creating a larger coupling zone.
Analysis of the present-day force balance (Lamb, 2006;Lamb & Davis, 2003) and geodynamic modeling of the long-term evolution of active continental margin in the Andes (Sobolev et al., 2006;Sobolev & Babeyko, 2005) and global subduction zones (Sobolev & Brown, 2019) suggest that sediments in subduction channels lubricate the subduction interface by decreasing the static friction coefficient. Sobolev et al. (2006) estimated the relation between the thickness of sediments at trench and the static friction coefficient of the subduction channel. To do so, they combined results of their geodynamic models with geodetic data at the locking depth of the subduction interface in Central and South Chile (Khazaradze & Klotz, 2003) and data on the thickness of sediments in the trench (Hoffmann-Rothe et al., 2006). They infer that the effective static friction coefficient is very low (about 0.01-0.015) if the thickness of the sediments in the trench is higher than~1.5 km increasing to about 0.05 or higher when the sedimentary layer thins to less than 0.5 km. Moreover, estimations of the static friction coefficient in the subduction fault as attested by observations of the heat flow by Gao and Wang (2014) propose friction in the range of 0 to 0.13, and smaller than 0.05 for subductions hosting GEQs. Using a similar approach, Hass and Harris (2016) calculate friction for Costa Rica subduction to be in the range of 0 to 0.12, with a preferred value of 0.06.
Since the suggested effect of sedimentary fill is to smooth the surface of the fault, it makes sense to also consider the incoming seafloor roughness. The studies of Lallemand et al. (2018) andvan Rijsingen et al. (2018) that consider seafloor data in the vicinity of the trenches of all subduction zones around the Ring of Fire show a statistical correlation between areas of smooth seafloor and the occurrence of great and GEQs. Observations demonstrate that low roughness of the seafloor at long wavelengths (80-100 km) is a crucial criterion for a subduction zone to host GEQs (supporting information Figure S1). The statistical correlation of GEQs with a smooth seafloor at long wavelengths might suggest that the roughness of the seafloor can limit the propagation of rupture (consuming energy of the rupture) and therefore limit the potential magnitude or even prohibit rupture nucleation at all. The subducted flank of ridges or sea mounts may prohibit propagation of the rupture along the trench axis, hence limiting the magnitude of the earthquake. This was supported by analogue modeling of earthquakes in the experiments of van Rijsingen et al. (2019) that modeled two types of surfaces in shear zone: a surface with small-scale roughness and a surface with large-scale roughness characterized by the occurrence of homogeneously distributed seamounts all over the length of asperity. The model with the rougher interface generated earthquakes with smaller rupture areas, lower magnitudes, and shorter recurrence times.
Another interesting observation is that GEQs occur mostly in subduction zones with neutral upper plate strain (in long-term), less frequently in the zones with compressive upper plate strain and never at zones of extensive upper plate strain (Heuret et al., 2012). The suggested interpretation of these observations was that moderate tectonic compressive stresses acting at the smooth plate interface should be high enough to allow frictional stresses to accumulate, but low enough so that rupture propagation is not inhibited (Heuret et al., 2012). Heuret et al. (2011) and Schellart and Rawlinson (2013) found weak correlations of the maximum observed Mw to the inclination angle of subduction and the width of the seismogenic zone. However, this correlation (pointed out earlier by Kelleher et al., 1974 andlater by Bletery et al., 2016) only becomes evident if events with magnitude higher than 8.5 are considered (Figures 1 and 2). In the last hundred years, earthquakes with magnitude exceeding 8.5 only occurred in subduction seismogenic zones wider than 75 km and dipping less than 35°; earthquakes with magnitude exceeding 9.2 only occurred in subduction seismogenic zones wider than 150 km and dipping less than 20°. All three mentioned subduction zone features, that is, dipping angles, sedimentary thicknesses (or long-wavelength roughness of the sea bottom), and strain regimes in the upper plate are combined in Figure 2. Note that most GEQs are hosted by subduction zones with low dipping angles, low roughness of the sea bottom, thick sediments in trenches, and with corresponding upper plates in a neutral strain regime.

GEQs and the Dipping Angle of Subduction Zones
Another parameter possibly affecting the maximum magnitude of earthquakes is the subduction rate (Ruff & Kanamori, 1980). Despite the fact that subduction rate did not show any correlation with the magnitude of the earthquakes in the analogue and simplified models of Corbi et al. (2017), it might still have an effect if we consider the temperature of the subduction interface. Advancing of the slab with a higher rate might cause the brittle-to-ductile transition to shift toward larger depths, which in turn will extend the width of the seismogenic zone. Thus, this aspect requires better investigation with numerical modeling.
The main problem in using the historical earthquake record to investigate parameters controlling GEQ magnitudes is that this record is too short, and in many regions incomplete. In other words, statistics of GEQs are limited. Therefore, to bridge this observational gap we turn to physics-based modeling. In this study we investigate the parameters controlling the occurrence and magnitude of GEQs, simultaneously trying to understand nature of the above-mentioned observed relations of the magnitudes of the GEQs with the dipping angles of subduction zones, smoothness and sediment cover of the incoming sea bottom, and with deformation style of the upper plates. To do so, we use numerical modeling of subduction earthquakes with a thermomechanical visco-elasto-plastic code, SLIM3D (Popov & Sobolev, 2008). The main advantage of our modeling approach is reproducing realistic stress and temperature distributions at the slab interface and in the overriding plate as well as using advanced, experimentally constrained nonlinear rheology in the subduction channel and elsewhere in the lithosphere-asthenosphere parts of the subduction system. We employ our recently developed cross-scale modeling technique (Sobolev & Muldashev, 2017) and design setups for subduction zones with various geometries, frictional parameters, and convergence rates. Our primary goal is to understand how the subduction angle, interface friction, and convergence rate influence the magnitude of earthquakes and deformation state of the overriding plate.

Method
Here we describe the principle setup and modeling approach of our method. In order to study the factors that control the magnitude of the largest possible events at subduction zones, we apply the cross-scale modeling technique from Sobolev and Muldashev (2017). This technique combines modeling of the subduction process for the millions of years time scale, as well as modeling of seismic cycles and earthquakes in 40 s to thousand years time scale using classical rate-and-state friction (RSF) rheology along the subduction interface. A similar approach using rate weakening rheology in order to mimic earthquakes at several years time scale was developed by van Dinther, Gerya, Dalguer, Corbi, et al. (2013) and.
We use the thermomechanical finite element numerical code SLIM3D (Popov & Sobolev, 2008), which solves the conservation equations of mass, momentum, and energy. The technique employs nonlinear rheology in which parameters for diffusion and dislocation creep are taken from laboratory rock experiments. The finite elements are isometric with a size of 3 km and our 2-D model domain is 900 km wide and 300 km deep.
First, we simulate the long-term subduction process with an advancing oceanic plate subducting under a continental plate over 5 to 10 × 10 6 years, with time steps spanning several thousands of years and dependent on the geometrical setup. We consider an oceanic subducting slab with a basaltic crust and harzburgitic lithospheric mantle, an overriding plate with a two-layer continental crust, a lithospheric continental mantle, and an asthenospheric mantle. For the shear zone of the subduction model, we consider the top layer of the elements of the subducted oceanic crust down to a model depth of 100 km. This shear zone or "subduction channel" is considered to have a small constant effective friction between 0.01 and 0.06 (Gao & Wang, 2014;Sobolev et al., 2006), as well as a weak "wet quartz" dislocation creep rheology from Ranalli (1997). For rheological parameters used in this study, see Sobolev and Muldashev (2017) and tables in appendix.
In order to obtain a realistic distribution of stresses and temperatures, we start each model by pushing the oceanic plate against the continental plate thus initiating subduction. Between the plates, we introduce a weak zone, inclined with a specific angle that influences the slab dipping angle. We further introduce velocity boundary conditions at the left side of the model that causes "mantle wind" in the asthenosphere also influencing the slab dipping angle. The bottom boundary and the asthenosphere part of the right side of the model remain open. After advancing of the slab to the bottom boundary of the model, we fix the velocity boundary condition at the bottom boundary of the model until the slab ceases to evolve (neither its geometry nor the stresses change anymore). After we achieve steady state behavior of the slab, we fix the velocity boundary conditions everywhere except for the surface. In this manner, we obtain a model with stable, adequate stress and temperature fields. The phases of preparation of the model are schematically illustrated in supporting information Figure S2.
Next, after having established a subduction model with realistic distributions of stresses, temperatures, and desired geometry of the slab (see Figure 1 in Sobolev & Muldashev, 2017, and Figure S2 in supporting information), we employ an adaptive time step procedure that allows us to model the interseismic loading of the shear zone with time steps of 5 years, and to approximate the processes of rupture nucleation and propagation with the time steps of 40 s. Postseismic relaxation is modeled with time steps ranging from 40 s to 5 years, depending on strain rates in the shear zone.

RSF
We employ the classic aging RSF law by Ruina (1983).
where τ is the shear stress, c is the cohesion, P f is the fluid pressure, σ n is the normal stress, and μ, μ * , and μ * st are friction, effective friction, and effective static friction, respectively. The a and b are parameters for the RSF law, and a * and b * are effective parameters for the RSF law. V is the velocity of slip on the fault, V st is the quasi-static (or reference) velocity (the velocity at which μ * ¼ μ * st ), L is the critical slip distance, and θ is the state variable. In continua approximation of slip in the subduction channel, we assume in 2 that , where _ ε p is square root of the second invariant of the plastic strain rate tensor and _ ε p st is the same at quasi-static conditions, taken as 10 -12.8 s −1 in our calculations (Sobolev & Muldashev, 2017).
Following (Sobolev & Muldashev, 2017), we explicitly integrate the state Equation 3, assuming constant slip velocity at every single time step and using RSF parameters that have been calibrated on the Valdivia earthquake in 1960, namely, This set of parameters simultaneously provides the correct characteristic stress drop, magnitude, and recurrence time of the Valdivia earthquake. For all our models we adopt the values of effective RSF parameters that we have estimated for the Valdivia earthquake but vary the static friction coefficient in the range from 0.01 to 0.06. Modeled stress drops of the earthquakes vary between 4 and 7 MPa in agreement with the observations for great and GEQs (Seno, 2014).
Similarly to the case of the Valdivia model from Sobolev and Muldashev (2017), we consider two "endmember" options for the depth dependence of RSF parameters: (1) depth-unlimited rate weakening friction with constant (a − b)* and b/(a − b) parameters starting from the depth of 100°C isotherm in the channel (Blanpied et al., 1995). In this way, the seismogenic region will be limited by the brittle-ductile transition.
(2) Rate weakening friction with constant (a − b)* ¼ −4 · 10 −4 and b/(a − b) ¼ 5 parameters starting from the depth of 100°C isotherm in the channel and down to the depth where temperature in the channel reaches 325°C (Blanpied et al., 1995), and rate strengthening for the deeper shear zone, with constant (a − b)* and b/(a − b) parameters of 4 · 10 −4 and −1, respectively (in accordance with Barbot et al., 2009).

Numerical Implementation of RSF
The SLIM3D technique employs the Newton-Raphson method to converge the solution during global iterations. Using the previous time step's solution as the initial guess allows finding the new solution in geological time scale models without any problems. However, during simulation of RSF the dependence of the shear stress on the current slip velocity (Equations 1 and 2) induces oscillations of the solution through the changing value of the yielding stress. In order to prevent oscillations and to achieve convergence, we implemented a procedure that guarantees that slipping elements of the fault are always on the yield at each global iteration. For that, at each global iteration, we solve the following nonlinear equation for _ ε p that corrects local visco-elastic trial stress in the channel: where τ is the square root of the second invariant of stress deviator approximating shear stress in subduction channel, c is the cohesion, μ is the friction, P is the dynamic pressure approximating normal stress in the subduction channel, _ ε p is the plastic strain rate (_ ε p ¼ V 2dx ; dx is the size of the element and V is the slip velocity), G eff is the effective visco-elastic shear modulus (Popov & Sobolev, 2008), and Δt is time step value. Close to mechanical instability, this equation has more than one root. We find extrema of the function using the dichotomy algorithm and always take the solution with the highest plastic strain rate.

From 2-D Seismic Moment to 3-D Moment Magnitude
All our high-resolution models are 2-D. Therefore, to compare our modeling results with observations, we must extrapolate our models to 3-D. To do so, we use an empirical scaling law for the rupture length and moment Magnitude, Mw (Strasser et al., 2010): where L is rupture length.
Taking into account that where M S,3D is three-dimensional seismic moment, M S,2D is two-dimensional moment, G is shear modulus, W is rupture width, and D is displacement during earthquake, we obtain the following relation between moment magnitude and two-dimensional seismic moment: In the following we use this particular relation to extrapolate results of our 2-D models to the nature. Note that we obtain a relation between 2-D and 3-D seismic moments rather than between rupture length and rupture width, which varies for different magnitudes. The robustness and limitations of this approach are discussed in section 4.5.

Model Configurations
In the first set of models, we use the same age for the lithospheric plates and the same convergence rate (7 cm/year) to exclude effects that are not related to geometry and static friction. We consider five different geometries, with dipping angles of 15°, 21°, 25°, 30°, and 45°, where the dipping angle is defined as the angle of inclination of the line connecting the ends of seismogenic zone: from 100°C in the channel to 325°C (Blanpied et al., 1995). The long-term models, labeled by their subduction inclination angle, are shown in Figure 3a. For each particular geometry, we consider models with effective static frictions of 0.01, 0.02, 0.03, 0.04, 0.05, and 0.06, and with both unlimited and limited velocity-weakening RSF parameters in the shear zone. We also consider two extra models with an inclination angle of 15°, with slow (3.5 cm/year) and fast (10 cm/year) subduction rates. For each model setup we simulate at least 50 earthquakes (Figure 3b). Since the first several earthquakes vary significantly in magnitude and recurrence time due to redistribution of stresses after activation of RSF, we exclude these first 10 earthquakes from our analysis (and from Figure 3b). The remaining 40 earthquakes exhibit quite low deviation in seismic moments and Geochemistry, Geophysics, Geosystems recurrence times, typically not exceeding 2% and at the most less than 5% from an average value (Figure 3b). All modeling results are presented in Table S1 in the supporting information.

Effects of Dipping Angle and Effective Static Friction
The maximum and average seismic moments for the earthquake cycles for all subduction geometries and static friction coefficients are shown in Figure 4a. It is clear that it is the slab geometry that has the largest effect on the seismic moment of earthquakes: The smaller the dipping angle of the slab, the larger the seismic events. Static friction has a smaller, but still significant, effect. The seismic moments are larger in subduction zones with lower effective friction coefficients. It is important to note that the results do not depend much on the type of RSF depth dependency tested in this study (Figure 4b).

Effect of Subduction Rate
In another set of models, we study the effect of subduction rate on the earthquake seismic moments. For this purpose, we run additional models with a 15°dipping angle of the slab and with different convergence rates, for different static friction coefficients in subduction channels. We test the subduction rates between 3.5 and 10 cm/year, which cover almost the complete range of subduction rates on Earth (e.g., Schellart & Rawlinson, 2013). Figure 4c shows the seismic moments for the modeled events in subduction zones with the same geometry (dipping angle 15°) but different static friction coefficients and subduction rates (3.5, 7, and 10 cm/year). It is clear that models with different convergence rates output very similar results. The reason is that, at relatively high velocities typical for subduction on Earth, the depth of critical isotherms (300-400°C) in subduction channels is not very different (see Figure S3 in the supporting information). Note, also, that the effect of different subduction rates is partially compensated by shear heating, which is higher for higher subduction velocities. It is clear that for a subduction velocity as low as 1 cm/year, the downdip depth of the seismogenic zone and its width will be significantly lower. However, such low subduction velocity is observed in collision zones rather than at typical subduction zones and, therefore, it is not considered in this study. Therefore, we conclude that the subduction rate does not affect magnitudes of GEQs, in agreement with the results of analog and simplified numerical modeling (Corbi et al., 2017). What is affected by the subduction rate is the period between the events, which scales with the inverse of subduction velocity, again in agreement with previous modeling (Corbi et al., 2017).
Therefore, our modeling suggests that in the hierarchy of effects on seismic moment, the most important is slab geometry and the second is the static friction coefficient. The slab velocity does not have an important effect for the typical subduction zones.

Main Role of the Seismogenic Zone Width
In Figure 5 we plot maximum magnitudes for all models versus the rupture width. The figure shows that for all models, irrespective of tested parameters, the maximum magnitudes of events are controlled by one feature, the rupture width. This conclusion fully supports the concept by Kelleher et al. (1974) with the addition of the concept by Ruff (1989), which suggests that the key factor for a high magnitude event to occur is a large and coherent contact zone between the downgoing slab and overriding plate that could be ruptured in a single event. Note that similar results were also obtained by Herrendörfer et al. (2015), Corbi et al. (2017), and Brizzi et al. (2020) using simplified numerical and analog modeling. This means that under the assumption that rupture length scales with rupture width, the models demonstrate that maximum magnitudes of the earthquakes are mostly controlled by the factors that increase rupture width. These factors are (1) a low slab dipping angle (the largest effect) and (2) a low friction coefficient in the subduction channel (smaller effect). Subduction velocity has a negligible effect.
The effect of slab geometry on rupture width is evident; that is, low-angle dipping slabs have wide seismogenic zones and hence can host wide ruptures.
The effect of friction in the subduction channel is less trivial and at first sight seems counterintuitive, suggesting wider ruptures and larger earthquakes in the subduction zones with lower friction in their channels. The nature of this effect is that the friction coefficient of the subduction channel affects rupture width by changing the depth of the brittle-ductile transition. If the friction coefficient is lower, then brittle (frictional) yield stress in the subduction channel is also lower and therefore brittle deformation regime in the channel extends deeper. Consequently, the rupture can propagate deeper and rupture width increases.
Our results on the dominant effect of rupture width are consistent with Kanamori's scaling of the seismic moment of earthquakes with the rupture area and stress drop (Kanamori & Allen, 1986). As observed in nature (and in our models), the stress drop of giant and great earthquakes does not vary much (Kanamori & Allen, 1986;Seno, 2014). In this case and if the rupture length scales with rupture width, then the seismic moment (and magnitude) of the earthquake depends almost solely on rupture width.
Our modeling shows that the seismogenic zone width, and hence the rupture width and maximum seismic moments of events, do not change considerably when subduction velocity is relatively high and varies only by a small percentage for the slabs subducting with velocities between 3.5 and 10 cm/year tested in our models. As most of the slabs on Earth are subducting with such high velocities, we do not expect any substantial effect of subduction velocity on magnitudes of earthquakes.
Because slab dip appears to be the most important aspect controlling rupture width and magnitude of the largest earthquakes, it is worth mentioning the major geodynamic processes increasing slab dip. These are high velocity of the overriding plate forcing fast trench rollback (e.g., Christensen, 2000), mantle flow tending to flatten the slab (e.g., Ficini et al., 2017;Hager & O'Connell, 1978), and a low negative buoyancy of the slab due to its young age or incorporated oceanic plateau in combination with fast overriding by the upper plate (e.g., van Hunen et al., 2000van Hunen et al., , 2002. Recently, modeling study by Brizzi et al. (2020) suggested that a decrease of the slab dip may also result from the increasing sediment thickness on the incoming plate. Note, however, that this is not a case in Central and South Andes where the dipping angles of the slabs are similar despite that the thickness of sediments is very different (Hoffmann-Rothe et al., 2006).

Does Shear Stress in the Channel Affect Magnitudes of Earthquakes?
The concept by Ruff and Kanamori (1980), based on the idea that the event magnitude is controlled by the strength of the mechanical coupling between the plates, is not confirmed by our models. We note, however, that the term "mechanical coupling" was used but was not defined explicitly in the paper by Ruff and Kanamori (1980). Here, we consider "mechanical coupling" as a shear stress in the subduction channel. Figure 6a shows that magnitudes of modeled events are dependent on the magnitude of the average shear stress at the subduction interface, and the observed relatively weak dependence is opposite to the common expectation that higher moment magnitudes should correspond to higher shear stresses. In fact, higher moment magnitudes correspond to lower shear stresses. The reason is that the lower frictional shear stress in the channel extends the brittle deformation regime to greater depths, which allows rupture to propagate deeper and thus increases magnitude of the earthquake. This confirms again that it is not the shear stress at the plates' contact zone that controls the maximum magnitude of an event, but rather its Figure 4. (a) Maximum (blue circles) and average (green circles) seismic moments for the events versus static friction in the subduction channel in models with dipping angles 15°, 21°, 25°, 30°, and 45°and with depth-unlimited rate weakening. Axis on the right side shows magnitudes of 3-D events calculated from 2-D seismic moments using the scaling described in section 2.3. (b) Average seismic moments for models with depth-unlimited rate weakening (green circles) and with rate strengthening deeper than the depth of the isotherm of 325°C (red circles) versus static friction in the subduction channel in the models with dipping angles 15°, 21°, 25°, 30°, and 45°. Axis on the right side shows magnitudes of 3-D events calculated from 2-D seismic moments using the scaling described in section 2.3. (c) Average seismic moments and magnitudes for the models with dipping angle 15°and different subduction velocities (orange circles denote 3.5 cm/year, green 7 cm/year, and magenta 10 cm/year) versus static friction in the subduction channels. Axis on the right side shows magnitudes of 3-D events calculated from 2-D seismic moments using the scaling described in section 2.3.
area. Another consequence of the widening of the seismogenic zone (i.e., deepening of its lower end) is increase of the average normal stress (Figure 6b). However, there is no causal relation between high normal stress and high magnitudes of events. Moreover, the highest average normal stress is modeled for the steepest subduction zone with relatively narrow, but deepest seismogenic zone (Figure 3, zone with 45°dipping angle), where the maximum magnitude of events is relatively low (Figure 6b).

Relation Between Deformation of the Upper Plate and Magnitudes of GEQs
In our cross-scale models we can also study the long-term deformation regimes of the overriding plates for subduction zones with different geometries and friction coefficients. For this purpose, we run models for each geometry and static friction coefficient for a period of 200 kyr without RSF (i.e., without earthquakes sequences) and observe motion of a point on the overriding plate at a distance of 200 km away from the trench. Landward movement of the observed point denotes compression; that is, tectonic shortening of the upper plate and seaward motion denotes its extension. Modeling suggests that pronounced extension of the upper plate (with velocity of the observation point more than 5 mm/year seaward) should be observed only for very steeply dipping slabs, in agreement with the observations by Lallemand et al. (2005); otherwise, the strain regime should be either neutral or compressive. Note that the largest events are modeled in the subduction zones with the neutral strain regime in the upper plates. The reason for the close to  neutral strain in the upper plates in most of our models is the low friction in the subduction channel that results in low mechanical coupling between the subducted and overlying plates.
Predictions of our models are consistent with the data (Heuret et al., 2011(Heuret et al., , 2012 that GEQs are observed only in subduction zones with either neutral strain in the upper plate (most frequently) or compressional (less frequently) and never in the case of an extensional upper plate (see large and small circles in Figure 7). Indeed, according to our models, pronounced extension is expected only for the steeply dipping subduction zones, which may host much smaller events than subduction zones with small dipping angles. Moreover, steeply dipping slabs are lacking sediments in the trench, which likely prohibits large ruptures due to strong heterogeneity of the topography of the subducting plates (Wang & Bilek, 2014).

Maximum Magnitudes in Models Versus Observations
The effects of the slab's dipping angle and the subduction channel static friction on maximum moment magnitudes are summarized in Figure 4a. For each particular slab geometry, the largest magnitudes are expected for the lowest friction coefficient, in our case 0.01. We consider these magnitudes as predicted maximum magnitudes of the subduction earthquakes for each specific slab dipping angle. We compare our model prediction with all historically observed GEQs in Figure 8a. If our predictions are correct, then the observed events at subduction zones for all dipping angles must have magnitudes lower or equal to our predicted values. As we see from Figure 8a that is indeed the case. Most of the magnitudes of the observed events plot significantly below the predictions. The closest to the predicted maximum earthquakes were the M9.5 Chile (Valdivia) earthquake of 1960 and the M9.0 Kamchatka event of 1952.
Our model specifies that, under the assumption of a subduction channel static friction coefficient of 0.01, all subduction zones with smooth interfaces (lubricated by sediments) and a dipping angle of less than 30°can host giant events with magnitude above 8.8, and with a dipping angle less than 25°, events with magnitudes above 9.2 ( Figure 8b). In producing the map of Figure 8b we used our predicted maximum magnitudes from Figure 8a and additionally assumed that events with Mw ≥ 9.2 (red line) can occur only in subduction zones with long wavelength roughness (RLW) less than 300 m over the length of trench segments not shorter than critical length of 400 km, and events with 9.2 > Mw ≥ 8.8 can occur only in subduction zones with RLW less than 500 m (see Figure S1 in the supporting information and Figure 2b) over the length of trench segments not shorter than critical length of 250 km. The critical lengths were defined as one half of the average lengths of ruptures with the Mw 8.8 and 9.2 (Strasser et al., 2010). The assumption about critical wavelength roughness is motivated by the empirical data suggesting that events with Mw ≥ 9.2 are observed only at trench segments with RLW < 300 m and events with 9.2 > Mw ≥ 8.8 at trench segments with RLW < 500 m (see Figure 2b). In the supporting information ( Figure S4), we present the prediction map produced without any assumptions about critical RLW values. We consider this map to be less reliable but we show it to illustrate the effect of ocean floor roughness. As we see from comparison of Figure 8b and supporting information Figure S4, the assumptions about roughness play an important role. For instance, without these assumptions, almost the entire active margin of South America is expected to host events with Mw above 9.2, while in the map built with such assumptions, there are large regions of the active margin of South America where that is not the case.
We note that the prediction map of Figure 8b and especially the map in supporting information Figure S4 are theoretical modeling-based predictions built using a number of assumptions and therefore should be considered with caution. For instance, our estimations of maximum magnitudes are almost certainly the upper-bound estimates because they assume the lowest friction coefficients in the subduction channels.

Geochemistry, Geophysics, Geosystems
Interestingly, our modeling-based predictions appear to be in a good agreement with the maximum magnitude predictions based on different kinds of statistical analyses of observations (Kagan & Jackson, 2013;Schäfer & Wenzel, 2019;Schellart & Rawlinson, 2013). For instance, regions where the giant events (Mw > 8.5) are most probable according to Schellart and Rawlinson (2013) (their subduction segments with ranking 5 and 6) correspond well with our regions of great events. Kagan and Jackson (2013) suggest that most of the subduction zones can host earthquakes with maximum magnitudes between 9 and 9.7, which agrees very well with our estimate of the maximum magnitude range of 8.8-9.6. A recent statistical study by Schäfer and Wenzel (2019), based on extrapolation of observed earthquake distributions to all major subduction zones using machine learning, also produce similar predictions (see map on their Figure 5). We find remarkable this consistency of the estimates of maximum magnitudes obtained by completely different methods.
Our models also predict that magnitudes of the events should be lower for the subduction zones with higher friction coefficients (i.e., smaller sediment thickness). This prediction is consistent with the conclusions by Gao and Wang (2014) who estimated friction coefficients in subduction zones using geothermal observations.

How Valid Are Our 2-D Model Predictions for 3-D Events?
The most important limitation of our modeling is its two-dimensionality. Using only 2-D models prohibits the study of how barriers at subduction interfaces, or the roughness of the interface in general, can influence an earthquake magnitude. We extrapolate our modeling results to the 3-D case using a statistical scaling law, which implies that rupture length scales with rupture width. However, that might not always be the case. If there were a barrier to the rupture propagation at a distance to the hypocenter closer than few rupture widths, or a large curvature of subduction zone segment (Schellart & Rawlinson, 2013), then this scaling Note that magnitudes of all observed giant earthquakes plot below the curve connecting model predictions for magnitudes of the greatest possible events for the subduction zones of specific dipping angles. Dashed lines indicate predicted maximum magnitudes assuming that the largest events release only one half or one third of the accumulated seismic moment and the rest is released in smaller events. Blue-gray zone indicates the domain where the maximum magnitude points would plot if more than one third of the accumulated seismic moment were released by largest events. (b) Location of the subduction zones with the predicted maximum events magnitudes (Mw) of 8.8 to 9.2 (thick orange lines) and of more than 9.2 (thick red lines) based on the dipping angle and long-wavelength roughness of seafloor. Dotted lines indicate subduction zones, where predicted maximum magnitudes are less reliable because of the lacking data on long wavelength roughness of seafloor. Circles show the location of Mw ≥ 8.5 subduction interface earthquakes (area scales with magnitude, color denotes upper plate strain regime). would fail. In this case, the event magnitude would be smaller than predicted by our 2-D model scaled to 3-D.
By definition, the rupture width and slip in earthquakes modeled in the 2-D plain strain approximation, assumes uniform deformation along the infinite (in fact very large) length along the trench. If in fact rupture length is not large enough, then our 2-D model-based predictions would fail. In order to estimate the critical length of the seismogenic zone for which our 2-D models are applicable, we conducted a series of 3-D experiments (see supporting information Figure S5). We performed 3-D and 2-D models with the same trench-perpendicular resolution of 6 km. Resolution along the trench in the 3-D models varied from 12.5 to 50 km. In our 3-D simulations we varied the length (L) of the seismogenic zone in order to find the critical ratio of length to width of seismogenic zone (W), L/W, at which the average 2-D moments (obtained from 3-D modeling) are close enough to the results of the 2-D models. It is important to note that 3-D models with different resolutions gave similar results. They suggest that a seismogenic zone with length greater than~2.5 times the width (Figure 9) can host earthquakes with similar seismic moments per length unit (3-D seismic moment divided by seismogenic zone length) as 2-D seismic moments predicted in 2-D models. As we see from Figure 9, the ratio between scaled to length 3-D seismic moment to 2-D seismic moment (M 3D /L/M 2D ) approaches about 0.9 at a L/W value of about 2.5, 0.95 at L/W of 3, and does not change significantly at larger values of L/W. Note that an overestimate of seismic moment by 10% leads to the overestimation of moment magnitude by only 0.03. According to Strasser et al. (2010), L/W is typically larger than 2.5 for the earthquakes with the magnitude above 8.5. This means that our 2-D modeling approach, with extrapolation to 3-D according to the scaling by Strasser et al. (2010), is reasonable for the GEQs with magnitudes above 8.5.

Concluding Remarks
We used a 2-D cross-scale numerical modeling technique to model giant subduction earthquakes and their seismic cycles to identify factors controlling the maximum magnitude of events. We extrapolated our 2-D modeling results to 3-D using empirical scaling relations between rupture lengths and moment magnitudes of great and GEQs. We tested effects of the dipping angle of subduction zones, the subduction channel friction coefficient (related to the thickness of subducted sediments), the subduction velocity, and the depth distribution of the RSF parameters. We also analyzed the modeling results to find out the relations between the maximum magnitude of events, the shear stress in the subduction channel, and the deformation style of the upper plate.
Modeling shows that it is the slab geometry that has the largest effect on the maximum magnitude of earthquakes. The smaller the slab's dipping angle, the larger the maximum seismic magnitude. The static friction coefficient in the subduction channel has a smaller, but still significant effect. The seismic moments are larger in subduction zones with lower effective friction coefficients (thicker sediments) in the channels. These results do not depend much on the type of depth dependency of the RSF (b − a)* parameter. The results are essentially similar for the models with depth-unlimited rate weakening and rate strengthening deeper than the 325°C isotherm.
The modeling also demonstrates that, under the assumption of rupture length scaling with rupture width, the maximum magnitudes of the earthquakes are exclusively controlled by the factors that increase rupture width. These factors are (1) a low slab dipping angle (the largest effect) and (2) a low friction coefficient in the subduction channel (a smaller effect).
The effect of slab geometry on rupture width is evident; that is, low-angle dipping slabs have wide seismogenic zones and hence can host wide ruptures. Figure 9. Ratio between scaled to length seismic moment (M3D/L) of the largest 3-D events to average 2-D seismic moments (M2D) for events in 2-D models versus the ratio of length to width (L/W) of the rupture zone. Blue dashed lines denote area of saturation of M3D/L/M2D, where this value becomes higher than 0.9. In all models, the subduction zone is dipping with an angle of 15°and the friction coefficient in the subduction channel is 0.015. In the 3-D models we select only events with maximum magnitudes rupturing the entire length of seismogenic zone. Blue solid circles denote models with trench-parallel resolution of 12.5 km, red circles ¼ 25 km, and yellow circles ¼ 50 km.
The friction coefficient of the subduction channel affects rupture width by changing the depth of the brittle-ductile transition. If the friction coefficient is lower, then brittle (frictional) yield stress in the subduction channel is also lower and therefore brittle deformation in the channel extends deeper. Consequently, the rupture can propagate deeper and rupture width increases.
Models suggest that there is a significant but not very pronounced negative correlation between maximum magnitudes of earthquakes and average shear stresses in the subduction channel. Our modeling also shows that subduction velocity in the range of typical subduction rates on Earth (between 3.5 and 10 cm/year) has a negligible effect on magnitudes of subduction earthquakes.
In agreement with observations, our models suggest that the largest earthquakes should occur in subduction zones with neutral (most frequently) or moderately compressive (less frequently) deformation regimes of the upper plate. This is a consequence of the low dipping angles and low static friction coefficients in subduction zones with the largest earthquakes, rather than a reason for the largest earthquakes.
Based on our modeling results we predict the maximum magnitude of the subduction earthquakes as a function of the slab's dipping angle. Magnitudes of all the largest historically observed events are lower than or close to our predictions derived assuming effective friction coefficient in subduction channels of 0.01. The closest to the predicted maximum earthquake magnitudes were the M9.5 Chile (Valdivia) earthquake of 1960 and the M9.0 Kamchatka event of 1952. An interesting prediction of our model is that all subduction zones with smooth interfaces (lubricated by sediments) and with dipping angles of less than 30°can host giant events with magnitudes above 8.8. For dipping angles of less than 25°, the magnitude can exceed 9.2 reaching value about 9.6 for the dipping angle of 15°. These results are in good agreement with the maximum magnitude predictions for subduction earthquakes based on statistical analyses of observations.
Our models also predict that magnitudes of the events should be lower for the subduction zones with higher friction coefficients (i.e., smaller thicknesses of sediments) in subduction channels in agreement with the conclusions based on estimated friction coefficients in subduction zones using geothermal observations.