## 2 EEPAS Forecasting Model  
Although inspired by the \(\Psi\) predictive scaling relations (Figure 2), the EEPAS model does not involve the identification of precursory seismicity for individual major earthquakes. It treats every earthquake as a potential precursor of future larger earthquakes to follow in the medium term [13]. Depending on the magnitude, this period can range from months to decades. The model has a background component and a time-varying component. The background component is a smoothed seismicity model, with the spatial distribution depending on the proximity to the location of past earthquakes. It is, in principle, time-invariant, but it is updated at the origin time of each contributing earthquake. The time-varying component, based on the \(\Psi\) predictive relations, is obtained by summing the contributions from all past earthquakes after a starting time \(t_{0}\) and exceeding an input magnitude threshold \(m_{0}\). The expected earthquake occurrence rate density is a function of the time, magnitude and location denoted by \(\lambda\). For times \(t>t_{0}\), magnitudes \(m\) exceeding a target threshold \(m_{c}\) and locations (\(x\),\(y\)) within a region of surveillance \(R\), the total rate density takes the following form:  
\[\lambda(t,m,x,y)=\mu\lambda_{0}(t,m,x,y)+\sum_{t_{i}\geq t_{0},m_{i}\geq m_{0 }}\eta(m_{i})\lambda_{i}(t,m,x,y) \tag{1}\]  
where \(\mu\) is an adjustable mixing parameter representing the proportion of the forecast contributed by the background model component; \(\lambda_{0}\) is the rate density of the background model; \(\eta\) is a normalizing function and \(t_{i}\) and \(m_{i}\) are the origin time and magnitude of the \(i\)th contributing earthquake, respectively. The contributing earthquakes come from a larger search region, which needs to be big enough to include all earthquakes that might affect the rate density within \(R\). The contribution from the \(i\)th earthquake to the rate density is given by  
\[\lambda_{i}(t,m,x,y)=w_{i}f(t|t_{i},m_{i})g(m|m_{i})h(x,y|x_{i},y_{i},m_{i}), \tag{2}\]  
in which \(w_{i}\) is a weighting factor and \(f,g\) and \(h\) are the densities of probability distributions which are based on the \(\Psi\) predictive scaling relations (Figure 2). These distributions depend on the magnitude \(m_{i}\) of the contributing earthquake. Following the notation in [20], the magnitude density \(g\) is a normal density of the following form:  
\[g(m|m_{i})=\frac{1}{\sigma_{M}\sqrt{2\pi}}\exp{\left[-\frac{1}{2}\left(\frac {m-a_{M}-b_{M}m_{i}}{\sigma_{M}}\right)^{2}\right]}, \tag{3}\]  
in which \(a_{M}\), \(b_{M}\) and \(\sigma_{M}\) are adjustable parameters. The time density \(f\) is a lognormal density of the following form:  
\[f(t|t_{i},m_{i})=\frac{H(t-t_{i})}{(t-t_{i})\sigma_{T}\ln(10)\sqrt{2\pi}}\exp {\left[-\frac{1}{2}\left(\frac{\log(t-t_{i})-a_{T}-b_{T}m_{i}}{\sigma_{T}} \right)^{2}\right]}, \tag{4}\]  
in which \(H\)(s) = 1 if s \(>\) 0 and is 0 otherwise and \(a_{T}\), \(b_{T}\) and \(\sigma_{T}\) are adjustable parameters. If all other parameters are fixed, the mean of the time distribution is proportional to \(10^{\sigma_{T}}\)Therefore, \(10^{\sigma_{T}}\) can be regarded as a temporal scaling factor. The location density \(h\) is a bivariate normal density of the following form:  
\[h(x,y|x_{i},y_{i},m_{i})=\frac{1}{2\pi\sigma_{A}^{2}10^{b_{A}m_{i}}}\exp\left[- \frac{\left(x-x_{i}\right)^{2}+\left(y-y_{i}\right)^{2}}{2\sigma_{A}^{2}10^{b_{A }m_{i}}}\right], \tag{5}\]  
in which \(\sigma_{A}\) and \(b_{A}\) are adjustable parameters. If all other parameters are fixed, the area occupied by the location distribution is proportional to \(\sigma_{A}^{2}\). Therefore, \(\sigma_{A}^{2}\) can be regarded as a spatial scaling factor. The adjustable parameters are fitted to maximize the log likelihood of the target earthquakes in the region of surveillance over a fitting period (\(t_{s}\), \(t_{f}\)) and a magnitude range (\(m_{c}\), \(m_{\max}\)).