## 3 Current Model and Numerical Algorithm
### Numerical Method
#### 3.2.2 Predictor Step  
The discrete equations in the _burning area_ in the predictor step respond to a semi-implicit Euler scheme:  
\[\frac{e^{n+1/2}-\bar{e}^{n}}{\Delta t}+\alpha u^{n+1/2}=r^{n}+q^{n}, \tag{41}\] \[e^{n+1/2}\in G(u^{n+1/2}),\] (42) \[\frac{c^{n+1/2}-c^{n}}{\Delta t}=-g(u^{n+1/2})c^{n+1}. \tag{43}\]  
where \(e^{n+1/2}\), \(u^{n+1/2}\), and \(c^{n+1/2}\) mean the predicted value of \(e^{n+1}\), \(u^{n+1}\), and \(c^{n+1}\), respectively. The basic idea is to treat the linear term implicitly. The heat contribution of non-local radiation \(r\) and fire-spotting \(q\) depends on the dimensionless temperature \(u\) and solid fuel mass \(c\), so the non-local radiation term \(r^{n}\) and fire-spotting term \(q^{n}\) are evaluated explicitly at time \(t^{n}\). Even so, Equations (41)-(43) continue to be nonlinear due to the multivalued operator \(G\). However, the solution to this problem can be reduced to explicit calculations.  
The multivalued operator in Equation (42) is maximal monotone, and hence, its resolvent \(J_{\mu}=(Id+\mu G)^{-1}\) is a well-defined univalued operator for any \(\mu>0\). Moreover, the Yosida approximation [29] of \(G\), \(G_{\mu}=\frac{Id-J_{\mu}}{\mu}\) is a Lipschitz operator, and the inclusion in Equation (42) is equivalent for all \(\mu>0\) to the following equation:  
\[e^{n+1/2}=G_{\mu}(u^{n+1/2}+\mu e^{n+1/2}), \tag{44}\]  
or  
\[u^{n+1/2}=J_{\mu}(u^{n+1/2}+\mu e^{n+1/2}). \tag{45}\]Now, rearranging Equation (41):  
\[u^{n+1/2}+\frac{1}{\alpha\Delta t}e^{n+1/2}=\frac{1}{\alpha\Delta t}g^{n}+\frac{1 }{\alpha}(r^{n}+q^{n}). \tag{46}\]  
and taking \(\mu=1/(\alpha\Delta\tau)\) in Equation (45), we obtain  
\[u^{n+1/2}=J_{1/\alpha\Delta t}\bigg{(}\frac{1}{\alpha\Delta t}g^{n}+\frac{1}{ \alpha}(r^{n}+q^{n})\bigg{)} \tag{47}\]  
which it is equivalent to solving  
\[(\alpha\Delta t\ Id+G)u^{n+1/2}\
i\bar{e}^{n}+\Delta t\ (r^{n}+q^{n}). \tag{48}\]  
Thus, denoting \(\mathbf{\propto}^{n}=\bar{e}^{n}+\Delta t\ (r^{n}+q^{n})\), the value of \(u^{n+1/2}\) is given by  
\[\frac{\mathbf{\propto}^{n}}{1+\alpha\Delta t} \text{if}\quad\mathbf{\propto}^{n}<(1+\alpha\Delta t)\ u_{\upsilon},\] \[u_{\upsilon} \text{if}\quad\mathbf{\propto}^{n}\in[(1+\alpha\Delta t)\ u_{ \upsilon},(1+\alpha\Delta t)\ u_{\upsilon}+\lambda_{\upsilon}],\] \[\frac{\mathbf{\propto}^{n}-\lambda_{\upsilon}}{1+\alpha\Delta t} \text{if}\quad\mathbf{\propto}^{n}\in[(1+\alpha\Delta t)\ u_{\upsilon}+ \lambda_{\upsilon},(1+\alpha\Delta t)\ u_{\upsilon}+\lambda_{\upsilon}],\] \[u_{p} \text{if}\quad\mathbf{\propto}^{n}\in[(1+\alpha\Delta t)\ u_{p}+ \lambda_{\upsilon},\infty]. \tag{49}\]  
Once \(u^{n+1/2}\) has been obtained, we calculate \(e^{n+1/2}\) and \(c^{n+1/2}\) explicitly from Equation (41) and (43), respectively,  
\[e^{n+1/2}=\bar{e}^{n}-\alpha\Delta tu^{n+1/2}+\Delta t\ (r^{n}+q^{n}), \tag{50}\] \[c^{n+1/2}=\frac{c^{n}}{1+\Delta g(u^{n+1/2})}. \tag{51}\]  
Notice that Equations (47), (50), and (51) can be solved simultaneously in all active nodes, and hence, parallel computation can be used to reduce the computational cost. In actual fact, the loop over all actives nodes to compute \(u^{n+1/2}\), \(e^{n+1/2}\), and \(c^{n+1/2}\) was parallelized using the API OpenMP [45].  
In the _burned area_, only Equation (41) needs to be solved. Since there is no fuel in this area and there are no changes in enthalpy or fire-spotting, Equation (41) simplifies to  
\[\frac{u^{n+1/2}-u^{n}}{\Delta t}+\alpha u^{n+1/2}=r^{n} \tag{52}\]  
which can be solved explicitly:  
\[u^{n+1/2}=\frac{u^{n}+\Delta t\,r^{n}}{1+\alpha\Delta t} \tag{53}\]