## 2 Governing Equations and Numerical Simulations  
In the flow through an urban environment, disregarding the thermal and buoyancy effects, the compressible effects are expected to be small due to the relatively low velocities. A requirement for accurate simulations is the use of a numerical discretization which allows for an adequate representation of the complete range of scales of turbulence. The incompressible Navier-Stokes equations are described in the following, where \(\mathbf{u}\) is the velocity vector, \(p\) the pressure, \(t\) the time and \(\
u\) the kinematic viscosity. \[ \begin{array}{l}\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot \
abla)\mathbf{u}=-\
abla p+\
u\
abla^{2}\mathbf{u}.\\ \
abla\cdot\mathbf{u}=0.\end{array}  \tag{2}\]  
The well-resolved LES was performed using the open-source numerical code Nek5000 developed by [PERSON] et al. (2010) at the Argonne National Laboratory and based on the spectral-element method (SEM) that is proposed by [PERSON] (1977). In this code, the Galerkin approximation is used for the spatial discretization and third-order discretization in time is used with explicit extrapolation and implicit backward differentiation for the nonlinear and viscous terms, respectively. The SEM combines the geometrical flexibility required to discretize the urban geometry with the high-order accuracy of spectral methods. As discussed below, the geometry of the buildings and streets is discretized using hexahedral elements and the scale disparity of the turbulent flows is resolved using the Gauss-Lobatto-Legendre (GLL) quadrature within each element. The turbulent flow through the simplified urban environment simulated here is shown in Figure 1, where the level of detail in the instantaneous flow structures can be observed. Note that the length and velocity scales are the building height \(H\) and the freestream velocity \(U_{\infty}\), respectively. The computational domain spans from \(x=-10\) to \(10\) in the streamwise direction. In the wall-normal and spanwise directions, the computational domain spans from \(y=0\) to \(2\) and from \(z=-3\) to \(3\), respectively. The inflow is a laminar boundary layer at \(Re_{\phi^{*}}=450\) (which is the Reynolds number based on displacement thickness) and a volume-force tripping is used to trigger the transition to turbulence as in Ref. (2017). We use the outflow condition developed by [PERSON] et al. (2019) and the top boundary has a zero-stress condition in the \(y\) direction together with a prescribed freestream velocity in \(x\) and zero velocity in \(z\). Finally, periodic boundary conditions are imposed in the spanwise direction. Figure 1: Instantaneous visualization of the flow under study showing vortical structures identified with the \(\lambda_{2}\) method (2017). The shown isosurface is \(-50\) (scaled with the freestream velocity and the obstacle height) and the structures are colored in terms of the streamwise velocity from (dark blue) \(-1\) to (dark red) \(1.8\). The well-resolved LES approach relies on a relaxation-term (RT) filter, which accounts for the contribution of the smallest unresolved turbulent scales through a dissipative force [41]. This force is given by \(-\chi H_{f}(u)\), then it is added to the right-hand side of the momentum equation in (2) and it is defined by a high-pass spectral filter \(H_{f}\) and a model parameter \(\chi\). In these simulations, \(H_{f}\) affects the two last modes and \(\chi\) is set to 200 based on \(H\) and \(U_{\infty}\). This technique was employed by [PERSON] et al. [42] to simulate the turbulent flow around a NACA4412 wing section up to a Reynolds number based on freestream velocity and chord length of \(Re_{c}\) = 1,000,000. Note that this method was thoroughly validated in the context of the spectral-element method by [PERSON] et al. [43] by comparing a fully-resolved direct numerical simulation of a turbulent wing at \(Re_{c}\) = 400,000 with a well-resolved LES of the same flow case, achieving excellent statistical agreement.