## 3 Model and Method
### Applied Numerical Methods and Simulation Parameters  
Deformation at the lithospheric scale can be modeled as a quasi-static thermomechanical coupling process. To simulate subduction, we solve the three conservation equations for momentum, mass, and heat [28]. The computations were performed using ASPECT version 2.2.0. The Stokes equation is expressed as follows:  
\[-\
abla\cdot[2\eta(\varepsilon(\boldsymbol{u}))]+\
abla p=\rho g \tag{1}\]  
In this equation, \(\eta\) represents the coefficient of viscosity, while \(\varepsilon(\boldsymbol{u})\) denotes the strain rate, defined as \(\varepsilon(\boldsymbol{u})=\frac{1}{2}\big{(}\
abla u+\
abla\boldsymbol{u}^{ T}\big{)}\), where \(u\) is the velocity and \(p\) is the pressure. The acceleration of gravity, \(g\), is 9.8 m/s\({}^{2}\). The density, \(\rho\), is approximated using the Boussinesq assumption for small temperature variations: \(\rho=\rho_{0}(1-\alpha(T-T_{0}))\), where \(T_{0}\) is the reference temperature, \(\alpha\) is the thermal expansion coefficient, and \(\rho_{0}\) is the reference density at temperature \(T=T_{0}\). The equations for mass conservation (Equation (2)) and heat conservation (Equation (3)) are as follows:  
\[\
abla\cdot\boldsymbol{u}=0 \tag{2}\]  
\[c\rho\frac{\partial T}{\partial t}+c\rho v\cdot\
abla T=K\
abla^{2}T-\rho SL+Q \tag{3}\]  
where \(c\) is the specific heat, \(\rho\) is the density, \(v\) is the velocity, \(T\) is the temperature, \(t\) is the time, \(K\) is the thermal conductivity, \(L\) is the latent heat of melting, \(S\) is the rate of melting, and \(Q\) is the radiogenic heat production rate. The parameters used in the simulations (Table 2) refer to previous research [29]. The depth-dependent mantle solidus and liquidus temperatures used in our calculations are discussed below.  
\begin{table}
\begin{tabular}{c c c c c c} \hline \hline \multirow{2}{*}{**Model Name**} & **Initial Overriding** & **Initial Subducting** & **Initial Dipping** & **Initial Subduction** & **Initial Overriding** \\  & **Plate Length (km)** & **Plate Length (km)** & **Angle** & **Rate \(V\) (mm/a)** & **Plate Movement** \\  & & & & & **Rate \(U\) (mm/a)** \\ \hline Case 1 & 100 & 200 & 45\({}^{\circ}\) & 30 & 60 \\ Case 2 & 100 & 200 & 45\({}^{\circ}\) & 30 & 40 \\ Case 3 & 100 & 200 & 45\({}^{\circ}\) & 20 & 50 \\ Case 4 & 100 & 200 & 45\({}^{\circ}\) & 40 & 50 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Parameters of the geometric models.  
A computational model was established with two dimensions: 300 km in width and 200 km in depth. The model incorporated a range of kinematic velocities for the Halmahera volcanic arc, where (\(U-V\)x > 0) indicated conditions favoring plate retreat. The boundary conditions for the models were established as follows: both the top and bottom sides were free to slide, while the left and right sides were designated as velocity boundaries. Different subduction rates \(V\) were applied in relation to the overriding plate movement rate \(U\). The boundary temperature was determined using the half-space cooling model as outlined by [30; 31; 32]:  
\[T=T_{0}+(T_{1}-T_{0})\frac{Z}{dZ} \tag{4}\]  
here \(Z\) is the thickness. A constant temperature of 5 \({}^{\circ}\)C is assumed to be maintained at the upper boundary, with the western boundary being adiabatic and the lower boundary set to 1350 \({}^{\circ}\)C.