## 5 Numerical Details  
The global discretization of the computational grid is accomplished using finite differences on a staggered Cartesian Arakawa-C grid [55]. The equations are spatially discretized on a fifth-order differential upwind scheme [56], while a third-order Runge-Kutta scheme is used on the temporal discretization [57]. The equations are implicitly filtered by the discretization of the computational grid and the subgrid-scale processes are calculated by the 1.5-order Deardorff numerical scheme [58]. Thus, it is assumed that the energy that is transported by the subgrid-scale vortices is proportional to the local gradients of the mean quantities [59; 60].  
Figure 6: GCI error bars estimated from the fine to medium grids for Scenario A: (**a**) for the normalized velocities \(\mathrm{u_{x}}/\mathrm{u_{so}}\) and (**b**) for the normalized concentration of the pollutant at the axial position of 0.5 H from the rear surface of the building.  
The convergence criterion maintained on each simulation is kept below the \(10^{-4}\) value for each computed variable based on the error. The time step is automatically adjusted for the CFL constant value of 0.9.  
Every simulation is carried out until flow stationarity has been achieved. Figure 7 presents the normalized velocity components with respect to the free-stream velocity (u\({}_{x}\)/u\({}_{\infty}\), u\({}_{y}\)/u\({}_{\infty}\), u\({}_{z}\)/u\({}_{\infty}\)), at the point with Cartesian coordinates of X: 10 H, Y: 5.5 H, and Z: 0.5 H and for the period from 600 to 1400 s for Scenario A.  
As it is shown in Figure 7, the flow seems to have attained stationarity conditions from the time point of 700 (s) onward. The attainment of stationarity conditions is confirmed by finding out that the mean value and higher order statistics are independent of the time of initiation of the measurements. If this is valid for the mean value and the autocorrelation function, the process is said to be weakly stationary. In this respect, the mean, the variance, and the autocovariance with a time delay of 2 (s) values are computed and given in tabulated form below (Table 1):  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Averaging Period (s)** & **Mean Value** & **Variance** & **Autocovariance** \\ \hline
600–800 & 1.0806 & 0.0473 & 0.3426 \\ \hline
800–1000 & 1.0851 & 0.0538 & 0.3738 \\ \hline
1000–1200 & 1.0879 & 0.0530 & 0.4316 \\ \hline
1200–1400 & 1.0889 & 0.0524 & 0.4443 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Mean Value, Variance and Autocovariance for Specific Time Periods.  
These computed values show that the flow field is at least weakly stationary.  
However, although the flow field has reached stationary conditions it is possible the concentration field has not reached it yet.  
Figure 8 presents the pollutant dispersion inside the computational domain for Scenarios A and C at the position with coordinates of X: 7.5 H, Y: 5.5 H, and Z: 0.5 H for the period of 200-1000 (s).  
The mean, the variance, and the autocovariance with time delay of 2 (s) values of the pollutant concentration are computed and given in tabulated form below (Table 2), which demonstrates that a passive scalar attains stationarity earlier in time than the velocity field: