## 2 Methods and Models  
A three-dimensional thermomechanical model originally developed from Stag3D [30] was employed in this study. The subducted oceanic plate geometry followed Slab2 [15]. The subduction direction in the central Aleutians changed from normal to right-lateral strike-slip faulting along the surface of the Queen Charlotte-Fairweather fault system, which underwent strike-slip-dominated deformation [4; 9]. The current convergent velocity and direction were obtained from MORVEL [31; 32]. The data were deposited into the TPDC repository, which has been granted a persistent identifier [[https://data.tpdc.ac.cn/en/disallow/8b266d22-fea7-4259-9a5f-8](https://data.tpdc.ac.cn/en/disallow/8b266d22-fea7-4259-9a5f-8) ac0 bd9e7869/]([https://data.tpdc.ac.cn/en/disallow/8b266d22-fea7-4259-9a5f-8](https://data.tpdc.ac.cn/en/disallow/8b266d22-fea7-4259-9a5f-8) ac0 bd9e7869/) (accessed on 4 December 2020).  
The model is prescribed with a size of 1600 \(\times\) 800 \(\times\) 400 km (width \(\times\) length \(\times\) depth) (Figure 2), and the node number is 80 \(\times\) 80 \(\times\) 100. The incoming oceanic plate is composed of a MORB layer at the top with a thickness of 7 km and is underlain by an ultramafic rock layer [33]. Moreover, the downgoing lithosphere geometry was kinematically prescribed without slab surface deformation, and the density and viscosity of each layer (e.g., upper crust, lower crust, subducted slab, mantle lithosphere, and asthenosphere for continent and ocean, etc.) remains constant (e.g., Tables in Supplementary Materials). The spatial distribution of ocean seafloor ages followed that from EarthByte [34]. Viscous decoupling was not included in the megathrust due to the uncertainty of frictional heating [35]. The calculation used the plate cooling model [36] with a temperature of 1100 \({}^{\circ}\)C at the base of the lithosphere and a mantle temperature of 1400 \({}^{\circ}\)C [36; 37]. The slab H\({}_{2}\)O content was calculated for P-T conditions computed from phase diagrams (e.g., [33; 38]). The thickness of the oceanic lithosphere was specified according to its age [39]. Because subduction in south-central Alaska has occurred since at least 160 Ma [40; 41], we prescribed the subduction time to be at least \(\geq\)20 Myr to ensure that the oceanic lithosphere reached a steady state with time, where the temperature variation was less than 0.1% with a lapse in time of \(\geq\)5 Myr.  
The composite model settings, including the model configuration, initial and boundary conditions, and physical parameters, are illustrated in the Supplemental Information.  
The modeling results are constrained by observations of the surface heat flow from the global heat flow database [42] and the estimates of the Curie depth [43] along parallel profiles across the arc (Figures S1 and S2). Our model follows the trench temperature boundary of the plate cooling model [36; 37]. We tested the resolution and found that the average temperature variance was \(<\)1% with a maximum temperature variance of \(<\)1.9% between 80 \(\times\) 80 \(\times\) 100 and 96 \(\times\) 96 \(\times\) 100 meshes. We performed sensitivity tests to investigate the robustness of our modeling results and varied the mantle viscosity from 1.0 \(\times\) 10\({}^{19}\) Pa s to 1.0 \(\times\) 10\({}^{21}\) Pa s and the mantle density from 3250 to 3350 kg/m\({}^{3}\). We present the benchmark model results as deviations from the reference models (\(\Delta\)T and \(\Delta\)H\({}_{2}\)O) and show these results at different depth levels within the oceanic slab. The tests show that mantle density variations (\(\pm\)50 kg/m\({}^{3}\)) induce small temperature variations of \(\sim\)10 \({}^{\circ}\)C at depth.