## 2 Materials and Methods  
In this study, our models were originally developed from stag3d [9]. An anelastic liquid approximation and the equations of conservation of mass, momentum, and energy [10] are used based on finite difference method. Based on the plate age and motion field mimicked for subducted curved plates [10; 11], we construct three-dimensional (3-D) thermomechanical models to estimate the thermal state-based distribution of wedge earthquakes in Alaska. Subducted oceanic plate geometry was incorporated into the model using Slab2 [12]. The intraslab subduction velocity field used MORVEL [13; 14] and interpolated the entire domain for the Pacific slab. The ocean seafloor ages of EarthByte [15] are used as an input of trenchward boundary conditions based on the plate cooling model [16; 17] (Figure 1). The thickness of the marine lithosphere was specified according to oceanic lithosphere age [18]. We prescribe the subduction time to be 20 Myr to ensure that the oceanic lithosphere thrusts to reach a steady thermal state. The modeling results are constrained by the observations of surface heat flow from the global heat flow database [19] and the estimate from Curry depth [20] along the parallel profiles across the arc. The thermal model is prescribed with a size of 1600 \(\times\) 800 \(\times\) 400 km (width \(\times\) length \(\times\) depth); the model grids are 80 \(\times\) 80 \(\times\) 100. The model includes the updated 3D geometry [12] and subduction velocity of the Pacific plate underlying the North American plate. The mechanical boundary conditions are as follows: the top boundary is rigid (dirichlet) and the bottom boundary is permeable (neumann) in the vertical direction; lateral boundaries are also permeable (neumann) except for the lateral boundaries, which have prescribed a time-dependent subduction velocity. Viscous decoupling is not included on the megathrust since the frictional heating is not robustly identified from the observation of surface heat flow. The evolving thermal structure of the subducting lithosphere is obtained from the calculated result of thermal modeling, which depends on the properties of the slab geometry, the age of the subduction plate, and the history of the subduction zone.  
Based on [PERSON] et al. [21] and [PERSON] et al. [22], we use a P-T-wt %-facies database to estimate the fluid upwelling sources on Alaskan megathrusts. The pressure (GPa) at every grid point is obtained by converting its depths (km) through preliminary reference earth model (PREM) parameters. Using the temperature and pressure provided by the numerical simulation, we estimate each rock facies (the mid-ocean ridge basalt and ultramafic rock)domain and the corresponding water content (wt %) at every grid. Through the interpolation method, we obtain the intraslab water content distribution (wt %) at various depths. The water content value difference between two neighboring points divided by interval distance reflecting slab dehydration (wt %/km) is obtained at each point.  
In this study, we focus on specific observations of wedge earthquakes and their underlying thermal and hydrous states in the Alaskan tectonic complex. To investigate the way the subduction regime controls or affects the distribution of continental wedge earthquakes, we focus on how the dip and slab geometry influence wedge quake distribution. The wedge earthquakes detected from 1 January 2000 to 31 December 2010 are employed in this study. The historical M \(>\) 7 earthquakes are also included.