## 2 Materials and Methods
### Numerical Solution Methodology  
The Finite Volume Method (FVM) is used to solve the system of compressible Navier-Stokes equations with temperature and pollutant concentration equations.  
We divide the calculation area into cells of size \(\Delta x\times\Delta y\times\Delta z\). The cell centers have coordinates:  
\[x_{i}=(i-0.5)\Delta x,\;i=1,2,\ldots,N_{x}\]  
\[y_{j}=(i-0.5)\Delta y,\;j=1,2,\ldots,N_{j}\]  
\[z_{k}=(i-0.5)\Delta z.\;k=1,2,\ldots,N_{z}\]  
We use an explicit scheme of the first order of accuracy ([PERSON]'s scheme):  
\[\frac{\phi^{n+1}-\phi^{n}}{\Delta t}=RHS(\phi^{n}),\]  
where \(\phi^{n}\) is any of the desired variables \((\rho,u,v,w,T,C)\), and \(RHS(\phi^{n})\) is the right-hand side of the corresponding equation.  
The continuity Equation (1) for discretization has the form  
\[\frac{\rho_{i,j,k}^{n+1}-\rho_{i,j,k}^{n}}{\Delta t}+\frac{(\rho u)_{i+1/2,j,k} ^{n}-(\rho u)_{i-1/2,j,k}^{n}}{\Delta x}+\frac{(\rho v)_{i,j+1/2,k}^{n}-(\rho v )_{i,j-1/2,k}^{n}}{\Delta y}+\frac{(\rho w)_{i,j,k+1/2}^{n}-(\rho w)_{i,j,k-1 /2}^{n}}{\Delta z}=0.\]  
The discretized form of the momentum Equation (2) for the component \(u\) is:  
\[(\rho u)_{i,j,k}^{n+1}=(\rho u)_{i,j,k}^{n}-\Delta t\left[\frac{ \left(\rho u^{2}\right)_{i+\frac{1}{2},j,k}^{n}-\left(\rho u^{2}\right)_{i- \frac{1}{2},j,k}^{n}}{\Delta x}+\frac{\left(\rho w\right)_{i,j+\frac{1}{2},k}^ {n}-\left(\rho w\right)_{i-\frac{1}{2},k}^{n}}{\Delta y}+\frac{\left(\rho w \right)_{i,j+\frac{1}{2},k}^{n}-\left(\rho w\right)_{i,j-\frac{1}{2}}^{n}}{ \Delta z}\right]\] \[-\Delta t\frac{\rho_{i+\frac{1}{2},j,k}^{n}-P_{i-\frac{1}{2},k}^ {n}}{\Delta x}+\Delta t[\text{diffusion terms}]_{jk}^{n}+\Delta t[\rho fx]_{jk}^{n}\]  
Similar equations can be written for the velocity components \(v\) and \(w\) along the \(y\) and \(z\) axes, respectively.  
The discretization of the energy (3) temperature (5) equations is performed similarly.  
The concentration Equation (6) is given in the following form:  
\[\frac{(\rho C)_{i,j,k}^{n+1}-(\rho C)_{i,j,k}^{n}}{\Delta t}+[\text{convective members}]=[\text{diffusion terms}]+S_{i,j,k}^{n}.\]  
For an explicit scheme, the Courant-Friedrichs-Lewy (\(CFL\)) condition must be met:  
\[\Delta t\leq min\bigg{(}\frac{\Delta x}{|u|+c},\frac{\Delta y}{|v|+c},\frac{ \Delta z}{|w|+c}\bigg{)}\text{-}CFL,\]  
where \(c\) is the speed of sound, and \(CFL\) is the Courant number (\(CFL<1\)).  
The model uses a structured grid with spatial resolution of 10 \(\times\) 10 \(\times\) 10 m\({}^{3}\) per cell in a domain of 1000 \(\times\) 1000 \(\times\) 500 m, resulting in a total of 5 million grid cells. Such a high spatial resolution allows detailed resolution of urban-scale features such as building-induced wind shear, urban heat islands, and pollutant trapping in street canyons.  
The temporal resolution is set to 1 s, ensuring compliance with the Courant-Friedrichs-Lewy (CFL) condition for numerical stability (CFL \(<\) 1). This fine temporal resolution permits accurate simulation of fast-changing atmospheric processes, including eddy formation and short-term wind speed fluctuations.  
Compared to coarser models (e.g., with cell sizes of 100 m or greater), our high-resolution scheme significantly enhances the ability to detect micro-scale airflow variations and pollutant stratification, which are critical for cities like Almaty with highly variable topography and meteorological inversion layers.  
The model's robustness was confirmed through sensitivity analysis of Peclet numbers, indicating dominant advection with balanced diffusion behavior. This configuration minimizes numerical dispersion and ensures physical fidelity of the computed fields.