Linked 3D modeling of megathrust earthquake-tsunami events: from subduction to tsunami run up

E. H. Madden, M. Bader, J. Behrens, Y. van Dinther, A.-A. Gabriel, L. Rannabauer, T. Ulrich, C. Uphoff, S. Vater, I. van Zelst 1 Department of Earth and Environmental Sciences, Ludwig-Maximilians-Universität München, Munich, Germany 2 Observatório Sismológico, Instituto de Geociências, Universidade de Brası́lia, Brası́lia, Brazil 3 Department of Informatics, Technical University of Munich, Garching, Germany 4 Numerical Methods in Geosciences, Department of Mathematics/CEN, Universität Hamburg, Hamburg, Germany 5 Seismology and Wave Physics, Institute of Geophysics, Department of Earth Sciences, ETH Zürich, Zürich, Switzerland 6 Department of Earth Sciences, Utrecht University, Utrecht, The Netherlands 7 Now at: Institute of Geophysics and Tectonics, School of Earth and Environment, University of Leeds, Leeds, LS2 9JT,

Computational advances now allow earthquake modeling to capture rupture dynamics on complex faults or fault systems on the scale of megathrust events (e.g. Murphy et al. 2016;Uphoff et al. 2017;Murphy et al. 2018;Ramos & Huang 2019). In addition, such models can include physical processes non-linearly coupled to the earthquake dynamics, such as off-fault plasticity (e.g., Andrews 2005;Dunham et al. 2011;Gabriel et al. 2013;Roten et al. 2014;Wollherr et al. 2018) or thermal pressurization of pore fluids (e.g., Bizzarri & Cocco 2006;Noda 2008;Gabriel et al. 2020). Such methods have realized realistic, observationally constrained simulations of several recent earthquakes, including the highly segmented crustal 2019 Ridgecrest, 2016Kaikoura, 2010Haiti, and 1992 Lan-approximations to the earthquake induced uplift as initial conditions (for a review see e.g. Behrens & Dias 2015).
To set tsunami model initial conditions, earthquake-induced seafloor displacements often are determined using an analytical solution for displacement due to a uniform rectangular dislocation within a homogeneous elastic half space (e.g. Okada 1985). A dislocation model may be taken from a finite fault model constrained by data inversion (e.g. Ji et al. 2002;Mai & Thingbaijam 2014;Bletery et al. 2016) as done by Allgeyer & Cummins (2014) and Jamelot et al. (2019), or designed to test certain source parameters as done by Geist & Yoshioka (1996). Stochastic models of seismogenic tsunami generation (e.g. Davies 2019) can either specify static slip on the fault following Andrews (1980) together with the idea that the final slip distribution after an earthquake rupture is self-affine, or use earthquake rupture models in the presence of stochastic stress (Geist & Oglesby 2014). In the context of an early warning system in Indonesia, complexity in space and time is achieved by including a grid of slip patches that together comprise a complex source model (Babeyko et al. 2010). Goda et al. (2014) highlight strong sensitivity of tsunami wave heights to site location and slip characteristics, and also to variations in dip, in stochastic random-field slip models for the 2011 Tohoku earthquake.
Several approaches also incorporate seismic waves into tsunami models.  present a 2-step tsunami modeling method: first, the spatial and temporal evolution of seafloor surface height is determined from a seismic wave simulation; second, this vertical velocity is used to source a 2D, nonlinear hydrodynamic tsunami model. As an earthquake source,  use a series of dislocations derived from dynamic rupture modeling of a potential future earthquake rupture in the Nankai Trough, Japan by Hok et al. (2011). Maeda et al. (2013) perform tsunami models of the 2011 Tohoku earthquake using tsunami-coupled equations of motion solved by the finite difference method (Maeda & Furumura 2011); these incorporate seismic waves and seafloor displacements generated from a 3D kinematic earthquake source. duction dynamics and seismic cycles to initialize the 3D dynamic earthquake rupture model (Fig. 1).
Such TECSEAS models, bridging the time scales of tectonic (TEC) and seismic cycle (SEAS, (Erickson et al. 2020)) models, provide insight into the role of rheology, temperature, subduction dynamics, fault geometry, loading and evolution, including spontaneously evolving splay faults (e.g., van Dinther et al. 2013bDinther et al. , 2014Sobolev & Muldashev 2017;Tong & Lavier 2018;Dal Zilio et al. 2019;D'Acquisto et al. 2020;Brizzi et al. 2020;Preuss et al. 2020). Earthquake rupture dynamics (including nucleation, propagation and arrest) are controlled by fault stress, strength and geometry and the surrounding material properties (e.g. Kame et al. 2003;Gabriel et al. 2012Gabriel et al. , 2013Galis et al. 2015;Bai & Ampuero 2017;van Zelst et al. 2019). Although these initial conditions may be informed by laboratory and regional observations (e.g. Aochi & Fukuyama 2002;Aagaard et al. 2004;Murphy et al. 2018;Ulrich et al. 2020), they remain difficult to constrain. This is particularly challenging in complex fault systems with lithological and geometric heterogeneities (e.g. Wendt et al. 2009), which impedes megathrust hazard assessment and mitigation. Constraints are particularly lacking in locations where observational data is sparse, either because earthquakes have not yet occurred or instrumentation is poor. Setting the earthquake model initial conditions from a subduction model provides much needed constraint on the earthquake model initial conditions. This approach also ensures self-consistency amongst initial conditions. For example, the assigned fault stress and strength are consistent with the fault geometry and material properties on and surrounding the fault. It also ensures self-consistency between those conditions and long-term subduction. For example, the assigned fault stress, strength and geometry are consistent with crustal, lithospheric and mantle deformation over geological time scales. In contrast to highly simplified earthquake models, subduction-initialized 2D dynamic rupture models produce earthquakes with multiple rupture styles, shallow slip accumulation and fault reactivation (van Zelst et al. 2019).
Here, we extend the approach that initializes a 2D dynamic earthquake rupture with a subduction and seismic cycle model (van Zelst et al. 2019) to a 3D dynamic earthquake rupture. The linked initial conditions include a curved, blind fault geometry, heterogeneous fault stresses and strength, and spatially variable material properties. This linkage requires consideration of the incompressibility and visco-elasto-plastic, plane-strain conditions of the subduction model versus the compressible, elastic conditions of the earthquake model. The resulting 3D dynamic rupture is linked with the tsunami model through the time-dependent seafloor displacements, following the same methods as in the first two examples.
For each application, we outline the methods, describe the resulting earthquake and tsunami, calculate the efficiency of each earthquake in generating that tsunami, and compare the tsunami results with those generated by a source incorporating the time-independent seafloor displacements from the 3D dynamic earthquake rupture -tsunami modeling 7 end of the earthquake model. We emphasize that these applications demonstrate the capabilities of the modeling framework; future, more involved and complex applications will certainly result in further knowledge gain.

METHODS
The earthquake and tsunami computational models utilized here are open-source, use discontinuous Galerkin schemes, and are facilitated by highly optimized parallel algorithms and software. The methods for one-way linking between these models and for setting earthquake model initial conditions from a subduction geodynamic and seismic cycling code are outlined in the following.
As a note on terminology; we use 'computational models' to describe the computer programs discretizing the model equations and implementing the numerical workflow, and 'physical models' to describe the structural setups, governing equations and selected input parameters. A 'scenario' refers to the results achieved by a computational model according to a specific physical model. We attempt to only use 'model' when the use of the term is unambiguous.

Earthquake modeling with SeisSol
SeisSol is the computational model used to simulate 3D dynamic earthquake rupture and seismic wave propagation (see App. A1). It solves the seismic wave equation in velocity-stress formulation using a Discontinuous Galerkin (DG) scheme with Arbitrary high-order DERivative (ADER) time stepping: ADER-DG (Dumbser & Käser 2006;Käser & Dumbser 2006). ADER is an explicit time-stepping method that achieves the same approximation order in space and time, but without requiring multiple stages for high discretisation order, as in e.g. Runge-Kutta schemes. The computational domain is discretized on a tetrahedral mesh, which simplifies automatic mesh generation for complicated geometries and facilitates static mesh adaptivity. Fast time to solution within SeisSol is enabled by recent hardware-aware computational optimizations targeting supercomputers with many-core CPUs (Breuer et al. 2014;Heinecke et al. 2014;Rettenberger et al. 2016) and an efficient local time-stepping algorithm (Breuer et al. 2016;Uphoff et al. 2017;Wolf et al. 2020). SeisSol has been validated against several community benchmarks (De La Puente et al. 2009;Pelties et al. 2012Pelties et al. , 2014Wollherr et al. 2018;Gabriel et al. 2020) following the SCEC/USGS Dynamic Rupture Code Verification exercises (Harris et al. 2009(Harris et al. , 2018.
SeisSol is specifically suited to solve for rupture propagation along complex, 3D fault geometries.
Its scalability enables large and long dynamic rupture models. For example, scenarios of the 2004 Sumatra-Andaman earthquake with up to 220 million element meshes and order 6 accuracy in space and time ran in 13.9 hours on the SuperMUC phase2 supercomputer (86,016 cores) at the Leibniz Supercomputing Centre in Garching, Germany (Uphoff et al. 2017 sam(oa) 2 -flash provides parallelization in shared (using OpenMP) and distributed (via MPI) memory . It scales up to thousands of compute cores, with problem sizes that exceed one billion grid cells with dynamic adaptive refinement and coarsening of cells .
Efficient adaptive mesh refinement is based on tree-structured triangular meshes (see App. A2). Fig. 1 shows schematically linkage between a subduction, earthquake and tsunami model. As we discuss here, this linkage requires consideration of assumptions inherent to each computational model, model dimensionality, and the time and space scales efficiently spanned by the computational models.

Tsunami initial conditions
Linkage from an earthquake model to a tsunami model requires several considerations. The recorded extent and the sampling rate of the domain and time frame of the earthquake model must be high enough to represent the required wavelength and frequency bands. For the transfer between models, the unstructured output from the earthquake model is bilinearly interpolated to an intermediate uniform Cartesian mesh. Also, we find that the earthquake model domain must be at least one fault length in each direction from the fault in order to minimize boundary effects. In terms of temporal resolution, 3D dynamic earthquake rupture -tsunami modeling 9 we find that a 1 Hz sampling rate of the earthquake displacement field is sufficient, as it is much smaller than the typical temporal scale of a tsunami waves.
To build the tsunami source, the time-dependent variations in the 3D vertical and horizontal coseismic seafloor displacements are mapped into vertical perturbations of the 2D initial bathymetry field of the tsunami model, similar to the approaches in ), Wendt et al. (2009) and Ulrich et al. (2019b. The time-dependent displacement fields are given by the 3D vector (∆x, ∆y, ∆z).
The east-west and north-south horizontal components, ∆x and ∆y, are incorporated into the tsunami source by the method proposed by Tanioka & Satake (1996), which combines the vertical component of the displacement vector with vertical changes induced by the interaction of horizontal displacement components and bathymetry gradients: where b = b(x, y) is the bathymetry (increasing in the upward direction). This contribution has been shown to be important both in general (e.g. Lotto et al. 2019) and for specific events, for example the 2011 Tohoku earthquake and tsunami (Murotani et al. 2015). However, for a physical model that includes flat bathymetry, the inclusion of horizontal displacement components has no effect on the tsunami source. ∆b is time-dependent, since ∆x, ∆y and ∆z are time-dependent. The tsunami is sourced by adding ∆b to the initial bathymetry and topography of the tsunami model. The tsunami computational model is capable of accurately modeling gravity wave propagation as well as non-linear effects due to advection and shoaling. The fastest resolved waves in the hydrostatic tsunami model are gravity waves, characterized by their linear wave speed of where H is the mean water depth and g the gravitational acceleration constant. For the applications presented in Sec.s 3 and 4, the average water depth is 2000 m, thus, the maximum wave propagation speed is approx. 140 m/s. Seismic surface waves in the earthquake rupture models are much faster, approaching approx. 2500 m/s. These seismic surface waves feature large transient amplitudes of up to 1 m, which is of the same order as the static uplift at the end of the earthquake models. In nature, fast seismic surface waves at the elastic-acoustic interface are converted into infra-sound or damped in the weakly compressible water column as the ocean response becomes non-hydrostatic at short wavelengths. However, seismic surface waves from an earthquake model may lead to spurious gravity waves in the tsunami shallow water approach. In addition, since a rupture model is not required to span the same spatio-temporal dimensions as the tsunami model, trailing seismic waves may show up as artifacts in final seafloor deformation. Thus, a space-time Fourier-transform based filter that removes unwanted signals from the tsunami source is applied.

Earthquake initial conditions
In the second example presented in Sec. 4, we initialize the earthquake model according to information from a subduction model that extends the long-term geodynamic model of Gerya & Yuen (2007) to seismic cycles using seismo-thermo-mechanical models (van Dinther et al. 2013a(van Dinther et al. ,b, 2014. This 2D subduction model solves for the conservation of mass, momentum and heat in an incompressible visco-elasto-plastic medium (Gerya & Yuen 2007, see App. C). After modeling millions of years of subduction, spontaneous frictional instabilities are simulated by reducing the time step combined with a strongly slip rate dependent friction law. Slip, stress and spatial characteristics of these frictional instabilities resemble those of natural earthquakes, albeit at very low slip rates (van Dinther et al. 2013a,b). We refer to these instabilities as 'slip events' to distinguish them from dynamic earthquake ruptures modeled with SeisSol, which capture frictional failure along a pre-existing fault and the accompanying seismic wave emissions.
To initialize the earthquake model using the subduction model, we port information from a single slip event following methods similar to those for initializing a 2D earthquake model by van Zelst et al. (2019), but here extended to a 3D earthquake model that is then linked to a tsunami model. We must pre-define a 3D fault geometry in the earthquake model and do so by extending the 2D fault that evolves during yielding in the subduction model in the third dimension. In addition, the material properties, stress state and friction coefficients from the 2D slip event are extended into the third dimension in the earthquake model. The computational subduction model is 2D and assumes plane-strain. The 3D computational earthquake model is not restricted to plane-strain conditions. We here extend the subduction model stress field into the third dimension for the earthquake model by setting the out-of-plane shear stresses to zero and the out-of-plane normal stress component to be a function of the two in-plane normal stresses and Poisson's ratio ν. Taking the in-plane coordinates from the subduction model as x and z and the out-of-plane coordinate as y, this is: σ yz = 0.
We here adopt ν = 0.5, which is the Poisson's ratio used in the subduction model. This value arises from the assumption of incompressible rock, a valid assumption for modeling the long-term deformation at convergent margins.
Because the assumption of incompressible rock is not appropriate for dynamic earthquake rupture modeling, Poisson's ratio must be reassigned to ν < 0.5 when transferring the material properties 3D dynamic earthquake rupture -tsunami modeling 11 from the subduction model to the earthquake model. Lame's parameter is then calculated from the assigned ν and the shear modulus from the subduction model.

Data management in massively parallel simulations
The earthquake and tsunami computational models and their linkage requires input at several stages.
For example, output from the subduction model is used to set the initial conditions for the earthquake model. Also, bathymetry data and seafloor displacements from the earthquake model are initial conditions in the tsunami model. To manage such data in massively parallel simulations, we use ASAGI (pArallel Server for Adaptive GeoInformation), an open-source library with a simple interface to access Cartesian material and geographic datasets (Rettenberger et al. 2016) (www.github.com/ TUM-I5/ASAGI). ASAGI organises Cartesian data sets as a collection of tiles. For dynamically adaptive simulations, the parallel partitions that are assigned to each compute node may grow or shrink,as the mesh resolution is adapted. ASAGI automatically replicates or migrates the corresponding data tiles across compute nodes, which greatly simplifies the computing access to material or geographic data at a specific location.

RUPTURE
Here, we demonstrate methods and highlight key results for linking a generic 3D megathrust rupture on a planar fault to tsunami generation and propagation in a basin with flat bathymetry and inundation on a linearly sloping beach. We compare tsunamis sourced by two earthquake scenarios that differ only by their near-surface fault strength, which controls the propagation of slip to the trench and results in one blind and one surface-breaching rupture.
We find that the surface-breaching rupture exhibits higher maximum and average fault slip and higher maximum and average vertical seafloor displacements. These differences are reflected in differing initial tsunami peak heights and wave shapes. However, the difference in peak height diminishes during tsunami propagation towards the beach and the inundation patterns are similar in both scenarios; inundation occurs along the same stretch of the beach and has the same run-up. The width of the inundated corridor inland from the coast and the timing of inundation do differ between scenarios, however, reflecting differences in timing in displacement of the water mass and in the magnitude of water displaced.
We can here study tsunami sensitivity to source time dependence, complementing heuristic findings (e.g. Davies 2019). By comparing these results to tsunamis sourced by the time-independent displacements, we find differences in the maximum run-up distances and in the temporal behavior of the tsunami.
where τ s is the shear traction on the fault, τ n is the effective normal traction, c is the on-fault cohesion, and µ s is the static friction coefficient. The right-hand-side in Eq. 3 is the static fault strength. τ n and τ s increase linearly with depth ( Fig. 2a,b) and we assign a uniform µ s of 0.275 (Fig. 2c). The effective normal traction includes the effect of pore fluid pressure, which we set to be depth-dependent and near-lithostatic (P f = ρ f gz, where ρ f = −3000 kg/m 3 , g is gravitational acceleration and z is depth, Fig. 2a).
The two earthquake scenarios that we compare here differ only by the static fault strength near the surface, controlled by c (Fig. 2d). In Scenario A, c = 0.4 MPa everywhere below 15 km, but gradually  In both earthquakes, the average rupture velocity resembles a typical megathrust tsunamigenic earthquake, but not a slower 'tsunami' earthquake Kanamori (1972). The mean values are similar in both scenarios, at 3.5 km/s for the blind rupture and 3.7 km/s for the surface-breaching rupture.
Rupture velocity remains subshear relative to the 4.3 km/s S-wave speed in the surrounding material along most of the fault during the blind rupture, but transitions locally to supershear speed up-dip from the nucleation location and along the upper part of the fault in the surface-breaching rupture.
For both earthquake scenarios, the maximum and minimum vertical seafloor displacements over the entire earthquake occur at t = 56 s. In the blind rupture, these are 2.6 m and −1.0 m, respectively.
For the surface-breaching rupture, the maximum is larger at 3.3 m, but the minimum is comparable at −1.1 m. The average vertical displacement at this time is 0.6 m for the blind rupture and twice as large, at 1.2 m, for the surface-breaching event. The displacements continue to change after this time in both scenarios until they reach constant values. For the blind rupture, this occurs at ∼ t = 80 s, when spatial maximum is 1.9 m, the minimum is −1.0 m and the average is 0.9 m. For the surfacebreaching rupture, this occurs later at ∼ t = 95 s, when the spatial maximum is 2.6 m and the minimum is −1.2 m. However, the average final displacement is 0.9 m, which matches that for the blind rupture. In both scenarios, surface waves continue to propagate until the predefined end of the model run, which is set to t = 120 s for the blind rupture and at t = 124 s for the surface-breaching rupture.

Tsunami propagation and inundation model and scenario
The tsunami physical model includes a flat seafloor and a linearly sloping beach (Fig. 4). Hence, we define the bathymetry by where x 0 = 200 km is the beach toe, where the slope begins. The initial sea surface is flat ('sea-atrest') and located at z = 2 km, which means that the coastline is located at x = 240 km. Above this level, the water depth is set to zero. We refer to the sea surface height (ssh) as the deviation from the reference height of 2 km. The model domain, which ranges from x = −400 to 400 km and from y = −400 to 400 km, is set to be larger than the horizontal extent of the earthquake physical model to minimize model boundary effects. The minimum mesh size is 12.2 m.
The bathymetry perturbation field ∆b (Eq. 1), accounting for the time-dependence of the simulated displacement field, is used to initiate the tsunami model. The unstructured output from the earthquake model is bilinearly interpolated to an intermediate uniform Cartesian mesh with a resolution of 1000 m, which is used for the transfer between models. We also apply a space-time Fourier filter to remove unwanted signals from the tsunami source that are present in the output from the earthquake model (Sec. 2.3.1). This filtering approach is further discussed in Sec. 5.1. Fig. 4 shows ∆b at t = 102 s.
Snapshots of the simulated tsunami wave-field are shown in Fig. 5. Also, key characteristics of the tsunami sourced by the blind rupture (Scenario A) and the surface-breaching rupture (Scenario B) are summarized in Table 2.
Maximum initial ssh resembles the maximum coseismic seafloor uplift of 2.6 m from the blind rupture and 3.3 m from the surface-breaching rupture. As shown in Fig. 6, the (filtered) source displacements of the blind rupture in Scenario A produce a smooth wave while those in Scenario B produces more abrupt initial displacements of the water column, as discussed in Sec. 5.1. The crosssection at y = 0 and t = 120 s in Fig. 6a shows that the wave peak is initially higher for the surfacebreaching rupture source. At t = 1200 s, just before first inundation in both scenarios, the heights of the peaks nearest the beach are more similar (Fig. 5b). For the blind rupture source, the tsunami reaches a maximum of ssh = 3.0 m at t = 1470 s, while for the surface-breaching rupture source, the tsunami reaches ssh = 3.4 m at t = 1480 s, a difference of 0.4 m. By t = 1600 s, the approximate time of maximum inundation, differences in peak heights have diminished, while tsunami wave distribution still differs (Fig. 5c).
We track the wave heights on land by comparing time series at 10 km inland from the coast near 3D dynamic earthquake rupture -tsunami modeling 15 x = 240 km along cross-sections at y = −150, 0 and 150 km (Fig.4d). The highest wave height occurs at y = 0 km in both Scenario A and Scenario B. However, the wave heights are asymmetric due to the uni-directional earthquake ruptures. Higher waves occur at y = 150 km, the part of the coast that is closer to locations of larger fault slip and uplift in both earthquake scenarios. Lower waves occur at y = −150 km, the part of the coast that is closer to the earthquake hypocenter. The height of the tsunami wave at the coast from the blind rupture is 0.8 m higher than the maximum wave height near the source, though this difference is only 0.1 m for the tsunami sourced by the surface-breaching rupture ( Table 2).
In general, wave peaks in Scenario B appear delayed relative to peaks in Scenario A, which is at least partially due to the fact that the location of highest seafloor displacement is farther away from the coast in Scenario B. This delay is more pronounced at y = 150 km and y = −150 km than at y = 0 km ( Fig.4d), where it is approximately 100 s. The average velocity from t = 1000 s to t = 1100 s of the wave peak from blind rupture source in Scenario A is 157 m/s, faster than that for the surfacebreaching rupture source in Scenario B of 142 m/s. This is in contrast to the earthquake rupture speed, which is on average faster in Scenario B and features supershear episodes.
The difference in propagation velocity does not affect the time of first inundation, which occurs at almost the same time in Scenario A (t = 1210 s) and Scenario B (t = 1220 s) ( Table 2). Inundation maps for both scenarios are shown in Fig. 7a and 7b. In both scenarios, the waves reach a maximum runup of 73 m at the center of the beach (near y = 0). Away from the center, the runup is lower. The inundation area is slightly asymmetric and skewed toward y = 150 km in both scenarios. Fig. 7c shows that, even though the first arrival occurs at approximately the same time, inundation in Scenario B is delayed near the coast and laterally along the coast by up to 100 s relative to Scenario A. This is probably due to the delay in the wave peak in Scenario B (relative to Scenario A) shown in

Comparison to tsunamis sourced using the time-independent displacements
We now compare the tsunamis from both scenarios when the time-independent filtered displacements from the end of the earthquake model are used in the sources instead of the time-dependent displacements from throughout the dynamic earthquake rupture scenario. The time-independent seafloor uplift has a maximum of 1.9 m for the blind rupture in Scenario A versus 2.6 m for the surface-breaching rupture in Scenario B. This is lower than the maximum uplift during the entire earthquake at t = 100 s of 2.6 m and 3.3 m in these scenarios, respectively (Table 1). However, the average time-independent displacements are equivalent at 0.9 m for both scenarios.
In general, the tsunamis from the time-independent sources inundate the coast earlier than those from the dynamically sourced models. To compare them, we equalize the times of first inundation: the time-independent source from the blind rupture in Scenario A is shifted by 40 s and the source from the surface-breaching rupture in Scenario B by 60 s.  later at the coast along y < 0 and earlier at the coast along y > 0. So, coastal inundation occurs over a shorter time. This is consistent with the source from a single time-step, rather than over the entire lifetime of the earthquake rupture. Farther inland, this trend is less pronounced. These differences in wave timing are visible in (Fig 8) as well, which shows that the tsunamis from time-independent sources lag behind those from the time-dependent sources at y = 0 km and precede the time-dependent sources at y = 150 km. The tsunami from the time-independent sources also over predict wave height at y = 0 km.

Summary
Using our virtual laboratory, we simulate tsunamis sourced by two megathrust earthquake scenarios that differ only by their near-trench fault strength and, as a result, slip behavior. This results in one blind and one surface-breaching rupture that differ in fault slip distribution and rupture kinematics.
The surface-breaching rupture exhibits 70% larger average fault slip and 40% larger peak fault slip.
Despite this, the surface-breaching rupture causes a maximum seafloor displacement that is only 35% higher than the blind rupture and the spatial averages of the seafloor displacements at the end of both earthquake scenarios are equal. The difference in maximum displacement is reflected in different initial tsunami peak heights, but this difference diminishes during tsunami propagation towards the beach and at the time of maximum inundation at the coast at y = 0 km, the peak wave height differential is only 0.4 m. The tsunamis in both scenarios flood asymmetric regions due to the uni-directional earthquake rupture and inundate the same 400 km of the coast and exhibit the same run-up of 73 m.
3D dynamic earthquake rupture -tsunami modeling 17 However, the tsunami generated by the surface-breaching rupture inundates a 15 % larger area overall due to a wider inundation corridor inland from the coast.
For both scenarios, using the time-independent displacements in place of the time-dependent displacements in the tsunami source results in later arrival at the coast, but faster coastal inundation.
For the surface-breaching rupture source, using the time-dependent displacements also over-predicts run-up.

TSUNAMI SOURCED BY A SUBDUCTION-INITIALIZED DYNAMIC RUPTURE
We now use output from a 2D subduction model (Sec. 2.3) to build the 3D earthquake model (Scenario C). In this way, earthquake initial conditions are assigned self-consistently and the tsunami source reflects the conditions developed over long-term subduction and seismic cycling. This adds complexity to the earthquake source relative to scenarios A and B in the previous section.
The earthquake initial conditions are heterogeneous and characterized by stark material contrasts, e.g., between sedimentary and basaltic regions, as well as smaller scale heterogeneous conditions reflected in fault geometry, fault strength, friction drop and stresses. Pore fluid pressure is elevated at a ratio of 0.95 to the lithostatic stress. Also, the critical slip weakening distance D c here varies with depth to resemble the friction drop measured during the geodynamic slip event.
We use the same methods as in Sec. 3.2 to link the earthquake model to the tsunami model. The fault does not intersect with the surface, so the rupture is blind, but it efficiently generates a tsunami. The sides of the earthquake physical model are 1600 km and it extends to 500 km depth (Fig. 9a).
The fault from the subduction slip event is a curve that extends for 320 km in the fault dip direction to a depth of 100 km. Fault dip gradually increases with depth, ranging from 2.3 • to 34 • and averaging 14.8 • . To pre-define the 3D fault in the earthquake model, fault locations are taken from the subduction slip event every 500 m in the x-direction (distance from trench) at the depth (z-direction) of the maximum strain rate over the entire slip event and smoothed with a moving average filter. To expand this line to a 3D fault surface, we assume that it is uniform in the third dimension, along fault strike, and extend it for 400 km (Fig. 9a). The fault does not intersect the surface, neither in the subduction nor earthquake models. We spatially discretize this structural model with elements of an edge length of 1 km on the fault, 5 km within a mesh refinement zone surrounding the fault, and a maximum element edge length of 20 km (App. B and Fig. A1).
The material properties and stress state are taken from the start of the subduction slip event In the earthquake model, fault failure occurs according the stress resolved on the fault as: This is equivalent to Eq. 3, but expressed in terms of the absolute normal traction, τ n , and the effective static coefficient of friction, µ s , at the beginning of the subduction slip event. We use the term effective here because µ s accounts for the pore fluid pressure, which has a ratio to the lithostatic stress of 0.95 along the fault during this slip event. τ n increases with depth, while τ s reaches a maximum near 43 km depth and minima near the surface and at large depths ( Fig. 10a,b). We assign the frictional cohesion, c, to equal the bulk cohesion, C, in the subduction slip event (Fig. 10d). We do, however, increase c in the earthquake model along the part of the fault that cuts through the sediments, as discussed at the end of this section. The Cartesian stress tensors are initialized in the earthquake model and SeisSol determines the tractions on the fault. Hence, it is critical that the fault and stresses in the earthquake model are in the same location relative to one another as in the subduction model (App. D).
After the onset of failure, the earthquake model fault weakens according to a linear slip-weakening friction law (Andrews 1976) and µ s changes with slip to an effective dynamic value, µ d . If the fault location is in a velocity strengthening region of the subduction model, we assign µ d to equal the maximum effective friction reached at that location during the entire subduction slip event, which may be locally larger than µ s . If the fault location is in a velocity weakening region, we assign µ d to equal the minimum value effective friction reached at that location during the entire subduction slip event. The 2D slip weakening distance, D c , is assigned such that the friction drop is equivalent to that 3D dynamic earthquake rupture -tsunami modeling 19 along the 1D fault in the subduction slip event (van Zelst et al. 2019). As shown in (Fig. 10e), such constrained D c varies with depth.
In a pre-modeling step taken before running the earthquake model, we check for points of static failure according to the earthquake model criterion (Eq. 5) along a cross-section through the 3D fault shown in Fig. 11. The failure criterion is met in three locations: within the shallow sediments, at one isolated point at 74.7 km depth, and in the region of 40-43 km depth.
In the subduction model, the slip event spontaneously initiates at 40-43 km depth. This is a location where the correct brittle failure mechanisms also are active. We predefine 3D a nucleating patch here in the earthquake model centered at x = 267 km, y = 0 km, z = −41.5 km and with a radius of 1.3 km. Inside the patch we assign µ s = 0.019, equal to the minimum value of µ s in the subduction model inside this nucleating region. In order to restrict nucleation laterally, we then set µ s = 0.025 in the regions outside of, but at the same depths as, the nucleation zone. This is the value of µ s above 40 km depth. We discuss linkage related to nucleation further in Sec. 5.2.
In the subduction model, nucleation of a slip event in the shallow sediments and at 74.7 km depth is inhibited. The shallow sediments are always at plastic failure, but velocity strengthening allows continuous creep through time without nucleation of brittle failure. However, localized exceeding of the failure criterion in the earthquake model would lead to rupture nucleation. Therefore, we prevent failure by assigning c = 5 MPa in the sediments above 25 km depth, which is the value of c in the deeper basalt.
At 74.7 km depth, high temperatures ensure deformation occurs predominantly through dislocation creep in the subduction model. Again, isolated exceeding of the failure criterion in the earthquake model leads to rupture nucleation here. To prevent this, we assign µ s = 0.02, which is the value of µ s nearby.
We make one minor, additional adjustment to the earthquake model. Near the material contrast at 27 km depth, µ s and µ d are anomalously large due to interpolation inaccuracies near stark material contrasts (Fig. 10c). We relax µ s and µ d at these points to the values of material below, such that µ s = 0.025 and µ d = 0.0097.

Earthquake rupture scenario
This subduction-initialized earthquake in Scenario C initiates spontaneously within the nucleation patch due to the locally high τ s and locally low µ s inherited from the 2D subduction slip event Waves continue to propagate after this time and surface waves and waves trapped in the sediments are still present at the end of the earthquake scenario at t = 241 s ( Fig. 12b and 12c). On average, the rupture velocity is 2.1 km/s, which is slower than the average velocities of the ruptures in Scenarios A and B, but still higher than a 'tsunami' earthquake Kanamori (1972).
The maximum vertical surface displacement over the entire earthquake is 28.1 m and occurs at t = 100 s (Fig. 13). At this time, the minimum vertical displacement is −5.6 m and the average is 3.6 m. After the earthquake simulation ends at t = 241 s (Fig.s 12b and 12c), the maximum and minimum displacements are 15.7 m and −6.7 m, respectively, and the average is 3.3 m. All characteristic displacements are larger than in Scenarios A and B (see Table 1).

Tsunami propagation and inundation
The tsunami physical model in Scenario C is the same as that used in Sec. 3.2, but its spatial dimensions are adjusted to the larger earthquake model. The beach slope begins at x 0 = 500 km, the coast is located at x = 540 km, and the size of the domain extends from x = −600 to 600 km and y = −600 to 600 km (Fig. 13). As in Sec. 3.2, the tsunami source, ∆b, incorporates the time-dependent 3D seafloor displacements from the earthquake scenario following Eq. 1. The unstructured output from the earthquake model is bilinearly interpolated to a structured mesh at a resolution of 1000 m. As previously done, we apply a space-time Fourier filter to ∆b which is discussed in Sec. 5.1.
The seafloor displacements are symmetric about the x-axis during the entire earthquake (Fig.s 12c and 13), as expected based on the symmetrically centered nucleation region. Fig. 14a-d shows the field of the sea surface height (ssh) for different moments in time. Initially, ssh reflects the vertical displacement magnitudes from the earthquake (Fig. 14a), then the tsunami develops a circular wave propagating away from the source (Fig. 14b) until it reaches the coast at t = 2050 s. At this time, the wave at y = 0 reaches ssh = 24 m, 0.4 m lower than the maximum height near the earthquake source, but much higher than in Scenarios A and B (Table 2). Fig. 14c shows the wave just after the time of first inundation and Fig. 14d shows the wave at the approximate time of maximum inundation 3D dynamic earthquake rupture -tsunami modeling 21 at t = 2400 s. Cross sections through the sea surface along y = 0 are shown in Fig. 15a-c. On average, the peak of the wave travels at a speed of 148 m/s, a propagation speed that falls between those for the tsunami in scenarios A and B in Sec. 3.

Comparison to tsunami from a time-independent source
We now model the tsunami using only the time-independent filtered displacements from the end of the earthquake model, instead of the time-dependent filtered displacements from throughout the dynamic earthquake rupture in Scenario C. The tsunami from this time-independent source arrives at the coast 200 s earlier than the tsunami from the time-dependent source, so we shift the simulation results by 200 s to synchronize the scenarios to the time when the coast is first inundated. The maximum timeindependent seafloor uplift is lower than the maximum uplift during the entire earthquake, though the mean displacement values for both cases are much closer (Table 1). than the spatial characteristics, as may be expected. Fig. 17a and Fig. 17b compare the wave profiles from the time-independent and time-dependent sources at t = 420 s. Along y = 0, the wave pattern is similar for both sources, but the timeindependent source produces a 1.1 m higher peak and this peak is advanced ahead of that in the original source ( Fig. 17a). At y = 150, the time-independent source produces a wave peak that is 1.3 m lower than for the time-dependent source, but the two wave peaks are in similar locations (Fig. 17b). Fig. 17c and Fig. 17d compare the waves for the static and dynamic cases in time at the coast (x = 540 km).
At y = 0, the time-independent source again over predicts the peak wave height (Fig. 17c), but the peak from the time-independent source is again slightly under-predicted at y = 150 (Fig. 17d).

Summary
In comparison to the two scenarios in Sec. 3, the average dynamic stress drop and rupture velocity for this subduction-initialized earthquake are lower, but fault slip and seafloor uplift are larger. The fault does not intersect with the surface, so the rupture is blind. We find that large and slow slip on a curved fault significantly impact tsunami characteristics, particularly maximum run-up, even in a fully buried and low stress drop event. Use of time-independent instead of time-dependent filtered displacements in the tsunami source over-predicts wave peaks at y = 0, but under-predicts peaks away from this location. This results in an under-estimate of the width of the inundation corridor everywhere except in the central region inland from the coast. Run-up is also over-predicted. Temporal differences include delayed arrival at the central coast, but advanced arrival along the more distant coast and at the locations farthest inland.

DISCUSSION
We here discuss selected aspects of the presented methods for linking subduction, dynamic earthquake rupture and tsunami models. We generalize on the use of time-independent versus time-dependent displacements as tsunami sources, compare the earthquakes and tsunamis in Scenarios A, B and C with real events, and contrast the tsunami-generating efficiency of the earthquakes in all three scenarios. We end with a look forward.

3D dynamic earthquake model -shallow water tsunami model linking
Use of a hydrostatic shallow water tsunami model in the linked modeling chain allows for evaluation of not only tsunami generation and propagation through open water, but also inundation at the coast.
In this workflow, we remove trailing seismic waves and specifically surface waves with a space-time Fourier-transform based filter. Fig. 18 shows the filtered vs. unfiltered source at t = 249 s, the end of the dynamic rupture model in the subduction-initialized earthquake (Scenario C). Note that the general uplift is kept unchanged, while the waves characterized by high ratios of frequency to wavelength (fast propagating waves/short wavelengths) surrounding the uplift area are effectively damped. We note that  use a low-pass filter, which does not completely eliminate the seismic waves from the tsunami signal near the source. Although we leave a full computational analysis to future work, we additionally find that the discretized shallow water solver also acts as a filter, since the fast wave creates a short impulse of momentum that decays immediately after the wave has passed, leaving the water column basically untouched.
The linked modeling approach permits study of tsunami sensitivity to source time dependence in 3D dynamic earthquake rupture -tsunami modeling 23 scenario cases, complementing heuristic findings (e.g. Davies 2019). We compare tsunami sourced by the time-dependent, filtered bathymetry perturbation with those sourced by the perturbation considering the time-independent displacements at the end of each earthquake scenario. We find, as may be expected, that temporal differences are larger than spatial differences.
In all three scenarios, use of a time-independent source results in more constant arrival times at the coast; the coast near the hypocenter is inundated later and the more distant coast is inundated earlier than for the tsunami with the time-dependent source. This behavior is consistent with a simulation Use of a time-independent source in Scenarios B and C (with surface-breaching and subductioninitialized earthquakes, respectively) over-predicts run-up. In Scenario C, the width of the inundation corridor is under-predicted at all distances from the coast. This is not the case for Scenarios A or B, however, where the width of the inundation corridor remains relatively unchanged. This effect in Scenario C may be related to over-prediction of the central wave peak (at y = 0) and under-prediction of the wave peaks away from here.
We see that the highest seafloor displacements over the entire duration of the earthquake do not control the tsunami heights during propagation to first order, but may control the width of the inundation corridor inland from the coast. This is seen in Scenarios A and B, which have seafloor displacements that differ by 35 %, but exhibit a decreasing difference in peak wave height during propagation.
The narrower inland inundation corridor for the blind rupture reflects its lower maximum seafloor displacements. This is also seen when comparing tsunamis from time-dependent and time-independent sources. In all scenarios, the highest seafloor displacements are consistently higher than the maximum displacements from the end of the earthquake that are used in the time-independent source (by up to an order of magnitude in case of Scenario C , Table 1). However, the peak wave heights from the time-dependent and time-independent sources are similar. This suggests that it is the more comparable average displacements that control tsunami wave heights, because they control the volume of water displaced.
In Scenario C, the inundation area of the tsunami from the time-independent source is narrower than that from the time-dependent source (Fig. 16c). This is also seen when comparing the tsunamis from the time-dependent sources from Scenario A versus Scenario B (Fig. 7c) and when comparing the tsunamis from the time-independent sources from Scenario A versus Scenario B (Fig. 7f). We attribute the narrower inundation corridors to the relatively low transient displacements during the associated earthquakes. However, this effect is not pronounced when comparing the tsunamis from the time-dependent versus time-independent sources in Scenario A (Fig. 7g) or from the time-dependent versus time-independent sources in Scenario B (Fig. 7h). These types of differences, shown here in generic models, may be challenging to distinguish from field data, for which regional and data-driven adjustment of the scenarios may be required. increasing static strength with depth (Fig. 2). However, what the initial conditions in the subductioninitialized earthquake in Sec. 4 reveal are shear traction and static friction depth profiles that vary with both depth and material, with the most obvious change from sediments to oceanic crust at approx. 28 km depth (Fig. 10). There is also variability in the static and dynamic friction coefficients with depth. This is a far greater (and continuous) range of values than are provided by laboratory measurements on a select number of samples. It also provides context for different results that may emerge from field studies of a single subduction zone.

Subduction
Because we assume that the fault, material properties and all linked on-fault conditions are constant in the third dimension not provided by the 2D subduction model, this approach still simplifies conditions relative to nature. Nonetheless, the heterogeneity directly influences the earthquake dynamics and therefore the tsunami source. For example, although the maximum fault strength is similar in all three scenarios presented here, the heterogeneous initial conditions for the subduction-initialized earthquake result in larger fault slip, but lower average dynamic stress drop and rupture velocity. This may correspond to natural megathrust behaviour. In addition, we observe that seismic waves traveling through complex materials around the fault (as opposed to a homogeneous material) influence earthquake dynamics, affecting rupture style and shallow slip accumulation.
These linking methods also provide avenues for investigating earthquake nucleation. We here run 3D dynamic earthquake rupture -tsunami modeling 25 the entire seismic cycling phase of the subduction model and then select one slip event from the series of quasi-periodic events. We then identify the beginning of this event as the time step when two adjacent points are at failure due to local exceeding of the failure criterion. In 2D coupling by van Zelst et al. (2019), these two points constitute the nucleation line on a 1D fault that lead to a dynamic earthquake rupture. Along a 3D fault, we must laterally restrict this location and do so by creating a 2D nucleation patch centered on these points at failure. Assigning a lower effective dynamic friction coefficient inside than outside of this patch, with both values taken directly from the subduction model, sustained dynamic earthquake rupture is initiated. Future work may investigate the sensitivity of nucleation to patch radius, patch shape or other characteristics. Also, future exploration of the relationship between slip initiation in the earthquake model and strain rate or slip rate in the subduction model may provide insight into the nature of slip nucleation itself.

Realistic scenarios?
Here we compare the linked scenarios against observed events, evaluating not only if a reasonable and/or realistic modeled tsunami is produced from a particular earthquake source, but also if that modeled earthquake is itself a reasonable/realistic event. We discuss the effects of slip-to-the-trench, and compare the tsunami-generating efficiency of the earthquakes in all three scenarios.

Blind versus surface-breaching ruptures (Scenarios A and B)
Prior to the 2011 Tohoku earthquake, the compliance and velocity strengthening behavior of shallow subduction zone materials was thought to prohibit fault slip near the trench. Most finite fault inversions restricted shallow slip, as for example slip inversions for the 2004 Sumatra-Andaman earthquake (Shearer & Burgmann 2010). In this case, the large Indian Ocean tsunami originally was attributed to These considerations motivated the two tsunami sources in Sec. 3. We expected that the blind and surface-breaching ruptures in scenarios A and B would produce distinctly different tsunamis and the surface-breaching rupture would inundate a larger area than the blind rupture (Fig. 7c). Inundation area is 15 % larger for the surface-breaching rupture. However, despite differences in the earthquake slip distributions and maximum slip values, the difference in peak wave height at the coast is small, and coastal inundation extent and run-up are the same in both scenarios. We note that the seafloor displacements in both scenarios differ less than the slip, partially explaining this contrast.
To better understand this, we calculate the efficiency of the Scenario A and Scenario B earthquakes in generating the resulting tsunamis as following Lotto et al. (2017a): where h max is the maximum wave height near the source (here, taken at t = 56 s) and < s > is mean slip. In Scenario A, A = 0.7. In Scenario B, B = 0.5. Thus, although Scenario B includes shallow slip and larger maximum seafloor displacements that contribute to a higher initial wave height, the source is less efficient at producing a tsunami than is the source in Scenario A. Such trade-off may explain why the differences between the two tsunami events are limited. We also note that the highly complex Scenario C earthquake matches the efficiency of Scenario A at C = 0.7 (see Sec. 5.3.2).
The relative loss of efficiency in Scenario B may be due to inefficiency of shallow slip in this surfacebreaching earthquake.
The presented 3D scenarios lead to tsunami efficiencies that are considerably higher than for those

Subduction-initialized earthquake and tsunami (Scenario C)
The subduction-initialized earthquake in Sec. 4 is a M w 9.0 event. The efficiency of this earthquake in generating the modeled tsunami is calculated using Eq. 6. With h max is taken near the source at t = 100 s, C = 0.7, which matches the efficiency of the blind earthquake rupture in Scenario A.

Looking forward
One main advantage of linked earthquake-tsunami modeling is that the model complexity can be readily increased or decreased, depending on the hypothesis being tested. For example, realistic representations of complex topography and bathymetry are permissible in both the earthquake and tsunami computational models, which may be critical not only for inundation modeling, but also for modeling tsunami genesis and propagation in deep water (e.g. Salaree & Okal 2020).
Similarly, we here restrict the off-fault constitutive behavior of the earthquake physical models to purely elastic and employ a linear slip weakening friction law on-fault. However, recent 3D earthquake-tsunami models of the 2004 Sumatra earthquake reveal the sensitive trade-off between shallow fault slip and off-fault elastoplastic deformation in controlling the tsunami height (Ulrich et al. 2020). Future linked modeling can help to elucidate how these aspects of earthquake dynamics influence tsunami behavior by readily switching to, for example, alternative friction laws such as rate-and-state, elastoplastic off-fault and wedge behavior, or thermal pressurization of pore fluids.
Tsunami generation from landslides is simulated well by established software (e.g. Clawpack, Mandli et al. 2016), but landslide sources are not incorporated into the presented modeling chain.
However, this chain can provide a way to test the importance of earthquake dynamics in the tsunami When initializing the earthquake model from a subduction model, we must honor the plane-strain conditions of the 2D subduction model while mapping the stress field into the 3D earthquake model.
This most likely does not reflect conditions in-situ. Ideally, one would incorporate the dynamic earthquake rupture into seismic cycling in the subduction model. Methodological advances may enable linking with a 3D subduction model and working toward this two-way coupling between earthquake dynamics and long term behavior.
Coupled feedback mechanisms beyond one-way linking from earthquake to tsunami also may be analysed in future work. For example, the two-way interaction of the water column and the solid Earth may be considered. As the feedback of pressure differences in the water column to the solid bottom is small in the models presented here (1 m of water wave height corresponds to approx. 0.1 bar or 10 kPa), the assumption of a rigid crust may not affect rupture evolution. However, assuming a rigid crust has been shown to influence the tsunami in the far-field and lead to earlier arrival times The scenarios presented in this paper are accompanied by detailed input files and model information (see App. A). This provides sufficient detail for other modellers to run all or parts of these scenarios in their linked model setup and compare their results to these. In future work, these applications may also be useful for community-wide comparison of dynamic earthquake-tsunami modeling approaches and alternative linking methods.

CONCLUSIONS
Surprises with devastating consequences in past earthquake-tsunami sequences motivate a better understanding of the physical connections between subduction, earthquake dynamics and tsunami from genesis to inundation. Modeling approaches bridging physical parameters and processes across these temporal and spatial scales are suited to help advance such research. We here use 3D dynamic rupture models as tsunami sources by building a virtual laboratory using open-source earthquake and tsunami computational models. We present three scenario applications to demonstrate the flexibility and capabilities of this linked modeling.
These methods are well-suited for hypothesis testing, such as isolating the influence of a single parameter on earthquake and tsunami behavior. In scenarios A and B, we compare tsunamis from a blind rupture and a surface-breaching rupture due to differing fault strength near the surface. The surface-breaching rupture does inundate a wider area inland from the coast, but we find that slip to the trench does not cause differences in inundation shape, change run-up or alter the length of coast impacted. This similarity may result from trade-off between the blind rupture's higher tsunamigenerating efficiency and the surface-breaching rupture's larger shallow slip, which leads to higher seafloor displacements and a larger displaced volume of water.
We use a subduction model to initialize the earthquake model in Scenario C. This approach provides reasonable earthquake initial conditions that typically are poorly constrained by data, but which exert first-order control over rupture behavior. Setting up the earthquake model in this way ensures physical consistency of the tsunami source with characteristics of both the subduction channel and the earthquake kinematics and dynamics. The resulting earthquake model includes a curved fault geometry and heterogeneous material properties and stress field, and the frictional parameters along the fault vary with depth. Due to these highly heterogeneous on-and off-fault conditions, the earthquake in this scenario has larger slip, but lower stress drop and slower rupture speed relative to the ruptures in scenarios A and B.
Future applications of these linked modeling methods can take advantage of the flexible adjustment of tsunami and earthquake model complexity. Direct studies of how subduction characteristics, earthquake initial conditions and earthquake dynamics govern tsunami behavior can help understand hazard in a given subduction zone. Such modeling may be specifically useful to constrain earthquake rupture and tsunami generation, propagation and inundation in complex megathrust systems, producing tsunami sources accounting for, for example, the effects of the slip to the trench, dynamic interaction between different fault segments (including splay faults), and off-fault coseismic deformation.

ACKNOWLEDGEMENTS
This effort has and continues to require intense cross-disciplinary collaboration and we want to acknowledge the strong team effort that has made it possible. In addition to the authors, a wide group has contributed to this effort. First, we would like to thank the Volkswagen Foundation (Volkswa-genStiftung) for extended funding and excellent support of the ASCETE and ASCETE II projects Wollherr.
Finally, participants in the 1st and 2nd ASCETE Workshops on Coupling Earthquakes and Tsunamis held in Sudelfeld Bayrischzell, Germany provided motivation for the development of these test cases and feedback on their design. We look forward to continuing discussions with the earthquake and tsunami modeling communities. We also would like to recognize the exceptional editorship by Editors Gabi Laske and Duncan Agnew, and to thank the reviewers, Joao Duarte, Brittany Erickson, Duncan Agnew and one anonymous reviewer, for their collegial and constructive reviews. We would like to thank the editorial team at GJI under Editor Prof. Gabi Laske for handling this cross-disciplinary manuscript in such a careful and constructive way. We also appreciate the collegial reviews from Joao Duarte, Brittany Erickson, Duncan Agnew and one anonymous reviewer.

DATA AVAILABILITY
We provide all necessary information and data sets required to reproduce the results presented in this work at https://tinyurl.com/yxn6zrqc. The archive includes all required configuration files, compilation parameters and input data. b inundation c propagation speed calculated for wave peak at y = 0 from t = 1000 to t = 1100 s, the time of first inundation       Zoom is to the region near the nucleation zone.   software at https://github.com/SeisSol/SeisSol. Documentation how to compile and use Seis-Sol is provided at https://seissol.readthedocs.io/.
SeisSol has been validated against a series of wave propagation and dynamics rupture benchmarks (see the references provided in Sec. 2.1) which are described in the "Cookbook" section of https: //seissol.readthedocs.io/.

A2 Tsunami modeling: sam(oa) 2 -flash
The sam(oa) 2 -flash tsunami model (Sec. 2.2) has been implemented within the sam(oa) 2 framework, which is publicly available as open source software at https://gitlab.lrz.de/samoa/samoa. Documentation on how to compile and use sam(oa) 2 and sam(oa) 2 -flash is provided at https:// samoa.readthedocs.io/. The sam(oa) 2 -flash version used in this paper is permanently archived in the gitlab repository, as branch samoa-flash-gji-2020. sam(oa) 2 -flash has been validated against a set of community benchmarks, following the test suite also used for validation of the underlying discontinuous Galerkin scheme, Stormflash2D (Vater et al. 2017). These benchmarks include tests for steady state solutions (Resting Lake) and the inundation of coasts (Oscillating Lake, Okushiri Tsunami). The archive at https://zenodo.org/record/ 3836668 includes instructions and required data files to reproduce these benchmarks. Results including a comparison to reference solutions are available at https://zenodo.org/record/3836668/ files/documentation_overview.pdf.
Tsunami modeling that includes inundation must handle varying spatial scales. sam(oa) 2 -flash uses adaptive mesh refinement (e.g., LeVeque et al. 2011), currently the most advanced approach for tsunami model meshing. This is in contrast to traditional tsunami modeling approaches employing a multitude of nested grids (e.g., Liu et al. 1995;Titov & Gonzalez 1997) or the alternative approach consisting of non-uniform, predefined meshes that use a priori knowledge about the wave behavior to assign high-resolution areas (Harig et al. 2008).

A3 Reproducibility
We provide all necessary information and data sets required to reproduce the results presented in this work at https://tinyurl.com/yxn6zrqc. The archive includes all required configuration files, compilation parameters and input data. such that the mesh generation may be run in parallel on a compute cluster. The scenarios presented here use mesh sizes of 13 million elements (Scenario C) and 16 million elements (scenarios A and B), which require computational resources well within the scope of typical applications for supercomputing centers or university clusters. Spatial resolution of earthquake faults has to be small enough to adequately resolve the dynamics behind the earthquake rupture front, where shear stress and slip rate vary significantly. The so-called cohesive zone, by analogy to tensile fracture mechanics, captures shear stress breaking down from its static to its dynamic value (Day et al. 2005). In the dynamic earthquake models used in scenarios A and B (Sec. 3.1), 400 m element edge lengths on the fault are combined with polynomial degree p = 5 (spatio-temporal order of accuracy of 6) leading to an effective numerical discretization distinctively higher. SeisSol's underlying numerical scheme defines initial conditions such as shear and normal stress at two-dimensional quadrature points located inside each tetrahedral element face which For the subduction-initialized dynamic earthquake model in Scenario C (Sec 4.1.1), we discretise the structural model with an element edge length of 1 km on the fault, 5 km within a mesh refinement zone surrounding the fault, and a maximum element edge length of 20 km. This properly resolves a minimum cohesive zone width of 193 m, which is appropriate along most of the fault. In Scenario C, it also is crucial that the meshed 3D fault matches the locations of the 2D fault from the subduction slip event (see Sec. D). We achieve this by extruding the 2D fault from the subduction odel along strike in GOCAD (www.pdgm.com), with straight horizontal lines between points at the same depth and no smoothing.

APPENDIX C: GEODYNAMIC SEISMIC CYCLE MODEL
The physical subduction model that we use has an extent of 1500 km in the x-direction by 200 km in the z-direction. The subducting oceanic plate consists of 4 km thick sediments, a 2 km thick basaltic upper crust, a 5 km thick gabbroic lower oceanic crust, and a lithospheric mantle ( Fig. A2a; Table A1).
The continental plate consists of a sedimentary wedge, an upper and lower continental crust and a lithospheric mantle layer. A constant velocity of 7.5cm/yr is applied to a small box inside the subducting plate to initiate and sustain subduction. At the sides and the top of the model, free slip boundary conditions are applied that only allow for tangential velocities along the boundary, while perpendicular velocities are not permitted. There is an open boundary at the bottom of the model where material can flow out, while free slip is externally applied (Gorczyk et al. 2007). The sticky air approximation, common in geodynamic modelling where topography develops, mimics an internal free surface, as low density, low viscosity 'sticky air' material is decoupled from the underlying rocks (Schmeling et al. 2008;Crameri et al. 2012).
The subduction geometry, lithological properties, temperature, viscosity, stresses and strengths develop spontaneously over 3.6 million years through solving thermo-mechanical conservation equations with a time step of 1000 years (compare Fig.s A2b and A2c). At temperatures below 100 • C, materials are velocity strengthening and a transition between 100 • C and 150 • C leads to velocity weakening at higher temperatures. The down-dip limit of the seismogenic zone develops as viscous deformation becomes progressively more dominant at temperatures above 350 • C (van Dinther et al. 2013b).
After 3.6 million years, a sufficiently steady-state subduction geometry has developed, suitable for a seismic cycle. The time step is gradually reduced by manually predefined intervals to 5 years, after which the seismic cycle phase of the model begins and lasts for approximately 30,000 years.
3D dynamic earthquake rupture -tsunami modeling 61 The subduction geometry (Fig. A2c) shows the oceanic plate subducting with an average dip of 14.8 • above 95 km depth during this phase.
After a spin-up period at the start of the seismic cycle phase, tens of quasi-periodic slip events occur (see van Zelst et al. (2019) , Fig 3). We refer to these as 'slip events' to differentiate them from earthquakes modeled with SeisSol. These slip events represent sudden, spontaneous increases of plastic strain rates during localized plastic failure in a narrow shear band (the 'fault' in the subduction model), which leads to reversed velocities due to elastic rebound and a drop of elastic stresses that were built up by the convergence of the subducting plate towards the overriding plate. Brittle failure is simulated according to Drucker-Prager plasticity and initiates when the second invariant of the deviatoric stress tensor σ II (Fig. A2d) meets the yield strength: σ II = C + µ sc 1 − (P f /P ) P (C.1) Here, C is cohesion, P is the mean pressure, P f is the pore fluid pressure and µ is the friction coefficient. Note that in the subduction model, there is no differentiation between bulk cohesion of intact rock and on-fault cohesion, c, after failure. We define the onset of a slip event as the first time step at which two adjacent points are at plastic failure (van Zelst et al. 2019). At this stage maximum slip rates are around 5.5 · 10 −11 m/s, with maximum slip rates of 9.3 · 10 −9 m/s reached during the entire event.
We know that this localized failure develops into a slip event because we select it from the series of events after running the entire seismic cycling phase of the subduction model (see Sec. 5.2 for further discussion).
Slip behavior after failure is viscoplastic rate dependent. In the strongly slip rate-dependent friction formulation (van Dinther et al. 2013a), the friction coefficient µ sc drops non-linearly from the static coefficient µ sc s to the dynamic coefficient µ sc d with increasing slip rate V , according to: where V c is the characteristic velocity at which half of the friction drop occurs.
Slip events occur mainly in the model subduction channel and the accretionary wedge. The events typically nucleate near the downdip limit of the model seismogenic zone in the basalt, after which they progress into the shallow sediments.
We choose one representative slip event to initialize the earthquake model. This event is chosen late in the simulation time of the seismic cycling, to ensure that the change of time step has no lasting effect on the slip events. In the chosen slip event, slip initiates at x = 220 km (according to the axis in

CONDITIONS TO THE EARTHQUAKE MODEL
To accurately map subduction model parameters to the earthquake model, it is crucial that the faults in both models are of the same geometrical shape. This ensures that they experience the same stress field and host the same on-fault properties. To verify this, before running the earthquake model, we compare fault locations, fault dip, effective shear traction, and failure on the 2D subduction model fault and along a 2D slice at y = 0 through the earthquake model mesh.In the following we show that the two faults have the same shape and capture the same initial conditions. 3D dynamic earthquake rupture -tsunami modeling 63