## 2 Materials and Methods
### Generalized Additive Model (GAM) Implementation  
A GAM was fitted to model vehicle counts as a function of spatial coordinates (longitude and latitude), day of the month, day of the week, and the spatial basis functions. The choice of GAM allows for including non-linear relationships between the predictors and the response variable (vehicle counts). The GAM model utilized Gamma distribution with a log link function. The model can be expressed as in Equation (1):  
\[\log(V(x,y,t))=\beta_{0}+\sum\limits_{i=1}^{7}\beta_{i}z_{i}+\gamma w+s(x,y)+s(t) \tag{1}\]  
where \(V(x,y,t)\) represents the vehicle count at the shelter/facility located at longitude, latitude pair \((x,\ y)\) on day \(t\); \(\beta_{0}\) is the intercept; \(z_{1},\ z_{2},\ z_{3},\ldots,z_{7}\) are (0, 1) dummy variables indicating the day of the week; \(w\) is a (0, 1) dummy variable indicating whether the school closure is in effect on day \(t\); \(\beta_{i}z_{i}\) represents the effect of each day of the week; \(\gamma w\) represents the school closure effect; \(s(x,\ y)\) and \(s(t)\) are smoothing spline functions for spatial coordinates and temporal dimension (days), respectively [39]. The GAM model (1) captures the evolving dynamics of evacuation-related shelter demand over space and time using spatial and temporal smoothing splines \(s(x,y)\) and \(s(t)\), respectively. The spline functions further allow to interpolate or forecast non-linear trends over space and time [39]. The categorical variable \(\beta_{i}z_{i}\) accounts for day of the week cyclical effects (i.e., differences in mobility patterns between weekdays and weekends) and the categorical variable \(\gamma w\) accounts for the effect of school closure on shelter demand over time. By utilizing these components our model is able account for deterministic temporal dependencies. However, if stochastic temporal dependencies in the data are also deemed to be significant, they can be incorporated in the proposed GAM model by including autoregressive terms [39], in a relatively straightforward way. The Generalized Additive Model (GAM) used in this study incorporates smoothing spline functions \(s(x,\ y)\) and \(s(t)\) to account for spatial and temporal variations in vehicle counts. These smoothing functions allow for the modeling of nonlinear relationships by balancing the trade-off between model flexibility and smoothness [40]. For the smoothing spline \(s(x,\ y)\) for the spatial coordinates, the following objective function is minimized:  
\[\min\limits_{s}\sum\limits_{i=1}^{n}\left(y_{i}-s(x_{i},\ y_{i})\right)^{2}+ \lambda\int\left(s^{\prime\prime}(x,\ y)\right)^{2}dxdy \tag{2}\]  
where \(y_{i}=\log(V(x,y,t))\). The first term ensures a good fit to the observed data, and the second term imposes a penalty \(\lambda\) to discourage roughness in the function \(s()\). The roughness penalty is quantified by the integral of the squared second derivative \(s^{\prime\prime}(x,\ y)\), which measures the curvature of the function. For the temporal smoother \(s(t)\) a similar minimization problem is solved. Therefore, the solution of the optimization problems yields the penalty terms \(\lambda=\lambda_{LatLon}\) for the spatial smoother and \(\lambda=\lambda_{Day}\) for the temporal smoother. The parameter \(\lambda\) governs the trade-off between these two terms [40] as follows:  
Low \(\lambda\) (High Flexibility): When \(\lambda\) is small, the penalty for roughness is minimal, allowing the function g to follow the data closely. In this case, the fit can become overly flexible, leading to potential overfitting. For example, the function might interpolate all data points, resulting in high variance but low bias [40]. Intermediate \(\lambda\): By carefully tuning \(\lambda\), a balance is achieved between fitting the data well and ensuring the function remains smooth enough to generalize to unseen data. This balance reflects the bias-variance trade-off inherent in statistical modeling [40].