On the Experimental, Numerical and Data-Driven Methods to Study Urban Flows

: Understanding the ﬂow in urban environments is an increasingly relevant problem due to its signiﬁcant impact on air quality and thermal effects in cities worldwide. In this review we provide an overview of efforts based on experiments and simulations to gain insight into this complex physical phenomenon. We highlight the relevance of coherent structures in urban ﬂows, which are responsible for the pollutant-dispersion and thermal ﬁelds in the city. We also suggest a more widespread use of data-driven methods to characterize ﬂow structures as a way to further understand the dynamics of urban ﬂows, with the aim of tackling the important sustainability challenges associated with them. Artiﬁcial intelligence and urban ﬂows should be combined into a new research line, where classical data-driven tools and machine-learning algorithms can shed light on the physical mechanisms associated with urban pollution.


Motivation
Urban areas are key elements in our society: around 75% of the population currently lives in cities in the European Union (EU) and it is expected that by 2050, 7 out of 10 people will be urban dwellers worldwide [1]. Cities are responsible for a significant fraction of the total greenhouse emissions in the world (above 60% according to Ref. [2]), and they are expected to play a prominent role towards the achievement of the climate-change targets from the Paris Agreement [3]. Two of the main challenges from cities for the coming years are the levels of air pollution and the so-called urban heat island (UHI) effect. According to the European Environment Agency (EEA), around 90% of the urban population in the EU was exposed to air pollution levels exceeding the ones recommended by the World Health Organization (WHO) between 2014 and 2016, and those levels of air pollution were responsible for around 500,000 premature deaths in 2015 [4]. In fact, more recent estimates indicate that this number may be as high as 800,000, where 40-80% of the cases are due to cardiovascular events [5]. On the other hand, the UHI phenomenon (whereby the temperatures in cities are higher than in the surrounding rural areas [6]) was related to 70,000 deaths due to the heat wave in Europe during the summer of 2003 [7].
Although the EU has introduced the use of predictive models to support pollutant concentration measurements [8], currently available techniques are unable to provide the spatial and temporal accuracy required to reproduce the pollutant-dispersion patterns within urban environments [9]. Similarly, more robust remote-sensing methods need to be developed in order to accurately predict the heat fluxes relevant to UHI phenomena [10]. Therefore, there is a pressing need for improved prediction and assessment methods to tackle these challenges and enable urban sustainability in the near future. Furthermore, it is necessary to gain further insight into the physics responsible for the pollutant-concentration and thermal distributions within cities. In particular, the coherent structures present in This article is structured as follows: in Section 2 we provide a review of the fundamental mechanisms present in urban flows; in Section 3 we give an overview of experimental campaigns aimed at understanding urban flows; on the other hand, Section 4 discusses a selection of numerical simulations of the flow in urban environments; an assessment of the relevance of coherent structures in urban flows, together with suitable data-driven methods, is given in Section 5; an outlook of the field and a summary of the main conclusions are provided in Section 6; finally, a selection of relevant modal decompositions to analyse urban flows is given in the Appendix A.

Fundamentals of Urban Flows
Urban environments can be described in terms of different length scales [16]: from the larger regional and city scales (in the range from 10 to 200 km), to the smaller neighbourhood and street levels (between 100 m and 2 km). Length scales are very important in the study of urban flows, since they determine the scope of the study and the possible techniques to be used. Intuitively, it is easy to see that one would study differently a street, i.e., a small scale, and a city, i.e., a large scale. At a very fine scale, e.g., within the range of 1 km, specific model is usually applied. The regional scale is affected by the urban area, and the flow responds to mean synoptic meteorological patterns perturbed by the city, as discussed by Britter and Hanna [16]. They also argue that at the city scale the "variations of the flow are mostly averaged out". On the other hand, most of the aspects affecting pedestrian comfort, air quality, emissions from vehicles, etc. need to be assessed at the street scale. In particular, the region of the flow directly affected by the obstacles is the urban-canopy sublayer [17]. Oke [18] provided an empirical description of the flow regimes visible in this layer, where three flow regimes are found depending on the ratio between the street width and the building height [19]: for narrow streets, the flow above the canopy can barely penetrate down to the street (skimming flow), and only one vortex is observed between the buildings; progressively wider streets lead first to the wake interference flow and later to the isolated-roughness regime, which exhibit much more interaction with the flow above the roof level and significantly more complex physics. In Figure 2 we show an instantaneous flow visualization in a simplified urban environment consisting of four arrays of blocks equally spaced in the streamwise direction and three in the span. This visualisation is the result of a high-resolution large-eddy simulation (LES) using the high-order spectral code Nek5000 [20]. The simulation uses a spectral-element mesh where elements are refined in the near-obstacle region. A Gauss-Lobatto-Legendre (GLL) quadrature of eight points is considered within each of the elements to ensure that the resolution is sufficient. Note that Figure 2 will only serve as a reference for the domain shape and dimensions. Please refer to Section 5 for a discussion on vortical structures and their interaction with the free-stream flow. This 4 × 3 model was simulated through well-resolved LES, i.e., a large-eddy simulation with very fine resolution, close to that of DNS. See more details in Section 4 and Ref. [21]. This figure illustrates the complex interactions of the vortical structures among the various arrays, together with the frame of reference considered for this work, as well as the geometrical parameters defining the urban environment. These parameters determine the regime observed in the urban area [18], in particular the width-to-height b/h and the height-to-separation h/ ratios. As shown in Figure 3, while b/h has a small impact on the resulting regime (particularly for tall buildings), the most critical parameter is h/ [18]. Instantaneous vortical structures identified with the λ 2 method [22] are shown with an isosurface of −40 (scaled in terms of the freestream velocity U ∞ and h). The structures are colored by streamwise velocity, ranging from (dark blue) −1.2 to (dark red) 1.8. Light grey indicates the bottom wall, and dark grey the buildings. Adapted from Ref. [18], with permission of the publisher (Elsevier).
The instantaneous flow in the three regimes discussed by Oke [18] is shown in Figure 4, where we consider two wall-mounted obstacles with varying streamwise separation. In this figure, also obtained through well-resolved LES [21,23], the flow is represented by vortical structures identified with the λ 2 method [22]. This method is based on the analysis of the the velocity-gradient tensor: denoting S and R as the corresponding symmetric and antisymmetric parts of the tensor (also known as the strain-rate and rotation tensors, respectively), Jeong and Hussain [22] define the vortex core as the region where the equation S 2 + R 2 has two negative eigenvalues. Note that it is sufficient to sort the eigenvalues in increasing order and evaluate the sign of the second one, a fact that gives the λ 2 name to this technique. All the obstacles have b/h = 0.5, and the various cases exhibit different values of leading to various interactions between the wake of the first obstacle and the leeward side of the second one. In Figure 4, we present the flow visualisation of three LESs, which aim at assessing the flow regimes described by Oke [18]. In this way, three simulations are produced, changing only the distance between the obstacles. Once again, those simulations were carried out using Nek5000 with a GLL spectral-element mesh. Roughly 205,605 elements were used to discretise the domain, i.e. a total of 105 millions grid points [23]. As far as flow conditions are concerned, a friction Reynolds number of Re τ = 280 is found in the vicinity of the obstacle [23]. All the simulations were carried out in the Cray XC40 system "Beskow", located at the PDC Centre for High Performance Computing at KTH Royal Institute of Technology. In the top figure ( Figure 4) we can visualize the skimming-flow regime, obtained with h/ = 1, in which the flow above the obstacles can barely penetrate in the region between the two blocks due to the very prominent effect of the first wake. A larger separation, leading to h/ = 0.5, is observed in Figure 4 (middle), which illustrates the wake-interference regime. Although there is still some interaction between the wake of the first obstacle and the second block, it is clear that some momentum transport from the flow above the canopy will take place between both obstacles. Finally, the isolated-roughness regime is observed in Figure 4 (bottom), with h/ = 0.25. Here the interaction between the two obstacles is limited, and therefore there will be significant momentum exchange between this region and the flow above the canopy.  [22], with an isosurface of −40 (scaled in terms of U ∞ and h). The flow exhibits three regimes depending on the obstacle separation, according to the classification by Oke [18]: (top) skimming flow, (middle) wake interference and (bottom) isolated roughness. The structures are colored by streamwise velocity, ranging from (dark blue) −1.2 to (dark red) 1.8.
The study of urban-flow regimes was expanded in the work of Zajic et al. [24], the authors of which analysed flow behaviour in a central business district (CBD). This study extended the description of Oke [18] to account for three-dimensional effects. In fact, important effects are found in the vertical edges which are only visible when the spanwise direction is considered. In three dimensions, flow analysis is segregated into the streamwise and spanwise coordinates. In the streamwise direction, i.e., following the x coordinate as shown in Figure 2, the length-to-height ratio w b /h is the principal factor affecting flow behaviour. For small values of w b /h the shear layer is completely detached from the obstacle. On the contrary, as this ratio grows the flow tends to exhibit a greater attachment to the surface and edges of the obstacle. An equivalent analysis can be applied in the spanwise direction y, this time as a function of the width-to-height ratio b/h . As the obstacle height grows, the recirculation bubble behind the obstacle tends to be overwhelmed by the side separation layer, resulting in intense turbulent fluctuations which may destroy the separation bubble [24]. The separation of buildings in the spanwise direction s changes the behaviour of the wake as well as the side-wall shear-layer interaction [24]. As the obstacles become closer in the spanwise direction, the wake-side-wall interaction increases. Additionally, the interaction changes the flow response to the surrounding physical factors. For instance, for an intermediate spanwise separation, the flow appears to be independent of the streamwise separation , due to the fact that the adjacent obstacle in the streamwise direction induces a jet behaviour on the flow. This phenomenon supports the need of a three-dimensional approach in the study of urban flows. Zajic et al. [24] also discussed the effects of non-uniform obstacle heights on the flow. In fact, larger differences in height lead to increased vorticity. This phenomenon results in more prominent turbulent structures.
Other works propose alternative classifications of flow regimes within the urban environment. Britter and Hunt [25] defined flow regimes in terms of the degree of interaction between the vortices and the freestream flow. According to Penwarden and Wise [26], the separation-to-height ratio of the buildings in the streamwise direction, i.e., /h (following the nomenclature introduced in Figure 2), can be used to establish the following criteria: for /h < 1 a stable vortex is found [25]; as /h tends to 1, higher speeds are found close to the ground. Note that the aforementioned behaviour is found when dealing with a uniform height over the domain. If significant differences between obstacle heights are found, the speed near the lower wall is significantly lower for the same separation-to-height ratio [25]. This last observation is supported by the work of Isyumov and Davenport [27], which shows that for large differences in height and a unitary separation-to-height ratio, the velocity near the lower wall is significantly smaller than what was predicted by Penwarden and Wise [26]. This discrepancy shows the importance of using three-dimensional approaches, since the various flow directions tend to behave independently.
Although the literature discusses certain observations regarding the flow features in these scenarios, there is a big need for more detailed information on the dynamical interaction of the flow structures in the various regimes. Furthermore, non-uniform street distributions and building heights will lead to a significantly more complex flow behavior, eventually responsible for the resulting distributions of pollutant concentration and temperature in urban flows. Those limitations suggest that understanding the actual behaviour of urban flows may require the implementation of more sophisticated approaches such as numerical simulations, especially in the study of the vortical structures.
Another important aspect of the study of urban flows deals with the characterisation of the atmospheric boundary layer (ABL). This is the region of the atmosphere that is in contact with the surface of the Earth. Although the ABL effects are sometimes ignored in the study of urban flows, in more specific applications it is critical to have a good characterisation of this boundary layer. Hlevca and Degeratu [28], together with Counihan [29], described the basic characteristics of a neutral ABL, i.e., without a thermal gradient. As with the classic turbulent boundary layer (TBL), the ABL is formed by two layers with an additional overlapping layer that lies between those. In the lower layer, i.e., 20-30% above the surface, the mean velocity exhibits a logarithmic behaviour [28] dependent on the roughness length, friction velocity and von Kármán constant. Note that the length roughness varies with the surface of the terrain under consideration and thus allows particularizing the law to the studied environment. However, Hlevca and Degenaru [28] explained that in the vast majority of engineering cases ABLs are modelled by means of a power law, thus approximating the ABL as a single layer. The power law is particularized to the application case by means of an exponent accounting for the terrain roughness. Equation (1) shows the power law presented by Hlevca and Degenaru [28]. It relates the ratio of mean streamwise velocities at two different heights (U 1 /U 2 ) with the heights z 1 , z 2 and two parameters particular to the domain: where the exponent α is particularised for a given domain and d is a reference distance particular to the domain. Counihan [29] explained the limitations of ABL models observing that various anomalies are found to exist in the published literature. Those anomalies are found in the definition of the specific parameters in both the logarithmic and power laws. Blocken et al. [30] showed that CFD simulations involving wall models tend to exhibit significant problems in the characterisation of the ABL. In fact, wall models typically incorporate wall-function roughness modifications based on experimental data, which become problematic when it comes to the alignment of the roughness length and height. The authors also present a solution strategy to avoid the aforementioned problem. It mainly consists of alleviating the requirements in the wall function. Nevertheless, these models tend to miss important information on turbulent structures. This explains why a part of the literature is focused on the study of turbulent processes within urban environments rather than the development of ABL models.

Experimental Studies
The complex, three-dimensional flow around an obstacle was thoroughly investigated using flow visualization by Martinuzzi and Tropea [31]. Later, Cheng and Castro [32] conducted detailed measurements based on X-wire anemometry combined with laser-Doppler anemometry (LDA) of the flow in the roughness sublayer within a simplified urban environemnt. Their results confirmed the strong three-dimensional (3D) character of the dominant structures, as well as by their interaction with the inertial sublayer. Further insight into the characteristics of the structures in the roughness sublayer (compared with those developing on smooth walls) was provided by Castro et al. [33] with additional measurements. These conclusions were complemented with spatial-correlation data obtained from particle-image velocimetry (PIV) by Reynolds and Castro [34]. On the other hand, the measurements by MacDonald et al. [35] focused on the characterization of turbulence statistics and on the development of scaling laws for modelling purposes. Note that the considerations for adequate wind-tunnel scaling of the full-scale urban environment are discussed in detail by Wang et al. [36]. The effect of building geometry on the flow field and pressure distribution were assessed experimentally by Xing et al. [37].
A wide range of experimental studies have attempted to identify the most relevant flow features, in terms of interaction between obstacles or influence of the incoming flow direction [38,39]. Although in Section 5 we provide a more thorough assessment of flow structures in urban environments, here we highlight some important findings from the experimental perspective. For instance, the interaction between the incoming flow and the horseshoe vortex forming around an obstacle was documented experimentally by Nagib and Corke [40] through flow visualizations, and this phenomenon was later analysed in detail through numerical simulation [41]. The importance of such structure knowledge for urban-boundary-layer parametrization was highlighted by Rotach et al. [42], whereas an attempt at a more explicit urban model was reported by Coceal and Belcher [43]. In this context, more realistic urban geometries were considered in the experiments by Rafailidis [44], while a shear-stress parametrization in a wind-tunnel experiment of a realistic urban geometry was proposed by Kastner-Klein and Rotach [45].
Detailed experiments of urban street canyons have been conducted in water channels with the aim of shedding light onto the dominant dynamic processes [46,47]. Other wind-tunnel studies have focused on the relevant implications of scaling down relative to the full-scale models, and they have provided detailed assessments of the possible sources of discrepancy [48,49]. Takimoto et al. [50] went one step farther, by evaluating the implications of scaling down a full-scale model in the context of coherent structures in the flow. Another element present in large-scale models is the disparity of building height, an aspect which was assessed in detail in the work by Hagishima et al. [51]. Other elements present in full-scale models which are challenging to realize and measure in wind-tunnel experiments are the effects associated to temperature variation. Kovar-Panskus et al. [52,53] discuss the impact of these effects in two detailed experimental studies. Furthermore, the presence of trees may reduce pollutant concentration due to the fact that they disrupt the streams present in the city canyon, as thoroughly documented in a number of experimental studies [54][55][56], although, depending on the canopy and building-block geometries, it can also increase pollutant concentrations (see [57]). Other elements, such as obstacles and obstructions in the streets, may have an important impact on the resulting flow, as discussed by Gayev and Savory [58]. Moreover, the potential of using modal decomposition to understand the role of the dominant flow structures was discussed by Kellnerova et al. [59], who performed proper orthogonal decomposition (POD) to shed some light into the data obtained with their PIV experiment. The role of modal decompositions in understanding the relevant flow structures in the context of urban flows will be discussed in Section 5 and in the Appendix A. Next, we will discuss both open-and closed-environment experiments of urban flows. Fernando et al. [60] gathered the data obtained in the Mock Urban Setting Test (MUST) [61], an experiment designed to study the flow in an idealised neighbourhood. The experimental site was conceived to mimic a uniformly distributed neighbourhood but using idealised geometries so as to limit the complexity of the flow. The experimental site is composed by a 10 × 12 array of semi-trailer containers located in the Desert Test Centre of the U.S. Army Dugway Proving Ground. The experimental site was conceived taking into account the variation in wind direction. In fact, the site exhibits very frequently south-north winds which induce a skimming-flow regime over the test area. The instrumentation used in the MUST experiment varied with the study considered. Fernando et al. [60] used a tower with sonic anemometers as main measurement system. Additional anemometers were installed at the front of the obstacle array. The experiment was modelled in a wind tunnel in order to compare the obtained results. While the MUST experiment mainly focused on flow characterisation, other urban-flow studies deal with the thermal effects. As discussed in previous sections, thermal effects are critical in the energetic behaviour of an urban area. Thus, further experimentation is required to characterize such effects. Fernando et al. [60] gathered a series of experiments providing data on the urban heat island (UHI), which is characteristic of many cities in the world as a result of the nature of urban environments and anthropogenic activity. In fact, UHIs have a direct impact not only on meteorological phenomena but also on several different areas such as energy consumption, air pollution, agriculture etc. [62]. The authors observe that the UHI has a direct impact on the local and regional meteorology of the surrounding areas, raising the temperature due to the heat exchange between the urban area (with a higher temperature and the adjacent zones) typically rural, which are naturally at a lower temperature. In fact, the authors show that in several of the analysed areas both the maximum and minimum temperatures in urban and rural zones raise over time. The temperature rise shapes the meteorology of the area, having a severe impact on the sustainability of both the city and the adjacent surroundings. Furthermore, thermal effects modify the structure of the turbulence within the urban area, which has a direct effect on the dispersion of pollutants and their interaction with the flow that remains trapped within the urban canopy.

Open-Environment Experiments
The Phoenix Urban Heat Island Experiment 2008 (PHEX-08) [63] gathered extensive temperature measurements aimed at validating UHI models. The PHEX-08 experiment analysed the circulation in street canyons due to temperature variation and how this circulation competes with the forced, i.e., vortex-induced, circulation that results from the interaction of the canyon fluid layer with the overlying ambient flow. The relation between both physical mechanisms, i.e., the balance between buoyancy and inertial effects, is described by means of the Richardson number. Temperature measurements were obtained by means of infrared (IR) thermography, sampling both the buildings in the street canyon and the surrounding areas. Furthermore, data was also obtained from a 300-m elevation. In the same line of work, other authors dealt with the UHI focusing on the underlying mechanisms and the parametrisation of the physical processes. Sharam et al. [62] explored the sensitivity of the UHI to urban parametrisations within the context of the Chicago metropolitan area (CMA) and its surroundings. A significant part of the available literature focuses on the development of physical parametrisations, i.e., models, that can help to identify the drivers of the UHI. Experiments tend to be used as verification or tuning methods for the developed models. In this case, Sharma et al. [62] gathered near-surface meteorological data from nine Mesowest stations located in the CMA to be compared with the simulations run during the study [62].
Moreover, Fernando et al. [60] presented the Hermoso Park dispersion study (HPS) which was carried out in Arizona to analyse the well-being of the residents in a 5-km 2 area. In this case the analysis focused on air quality and in particular the physical phenomena that affect it. The experiment covers measurements of meteorology, turbulence, PM 10 and PM 2.5 (i.e., particulate matter with diameter smaller than 10 and 2.5 µm respectively) among others. However, Fernando et al. [60] cover exclusively the flow, turbulence and transport data. The measurement system mainly consists of a meteorological tower equipped with sonic anemometers to capture the aforementioned quantities. Moreover, balloons were launched at the same location to measure velocity, temperature and pressure as well as to gather PM data. Furthermore, additional meteorological stations were distributed over the site in order to gather additional measurements.
As far as the results are concerned, Fernando et al. [60] used the gathered measurements to present a set of field observations mainly in terms of wind profiles, mean velocities, turbulent kinetic energy, temperature (and connected UHI effects) and pollutant concentration. They compared the thermal and pollutant results with a CFD model in order to analyse the suitability of using predictive CFD models for the study of urban flows. Their conclusions focused on the methodological perspective of the issue, arguing that both field measurements and models shall be combined. In fact, their observations showed that urban flows typically exhibit several physical phenomena which critically interact with one another, e.g., pollutant dispersion and thermal phenomena. Hence, the use of simplified approaches, such as controlled-condition experiments or highly restrictive models, exhibit severe flaws in the identification of such interrelated phenomena and by extension in the study of urban flows. For instance, environmental models tend to consider coarse resolutions-e.g., air-quality studies (U.S. EPA) use a resolution of 1 km-which are typically inadequate due to the strong gradients present in the buildings and street canyons [60]. As far as pure experimental approaches are concerned, confined-environment experiments tend to have problems in matching the actual physics of urban flows, mainly due to the difficulties of reproducing urban environments in a laboratory and due to the complexity of reproducing certain mechanisms of turbulence due to confinement effects. For example, Gadilhe et al. [64] showed how the measurements of velocities in the recirculation and wake regions of the urban canopy were not reliable.
Vita et al. [65] present a complete study on the assessment of pedestrian distress in urban environments comparing full-scale experimental techniques with wind-tunnel approaches as well as numerical methodologies. They distributed eight sonic anemometers over the streets at a distance from the ground of 2 m. In addition to those probes, a reference anemometer is permanently working. This reference probe was placed on the roof of a 62-m-height tower in a 10 m mast with the purpose of reading the baseline conditions in the undisturbed zone. Both reference and testing probes were configured to record three-dimensional data. Note that the testing zone was selected to be highly gusted, precisely because the objective was to assess the different methodologies used in pedestrian distress studies. The full-scale testing was aimed at obtaining data on meanflow speed. In addition, the data was also used to help characterising the flow. The measurements were then compared with the results of wind-tunnel experiments and numerical simulations. Comparing the measurements with the numerical results, the authors assessed the suitability of the methods in the study of turbulent flows within the context of urban environments. Mean fields, i.e., mean-flow velocities happen to exhibit a good agreement in the vast majority of the domain considered, when comparing measurements with numerical simulations. However, some areas such as the recirculation areas were particularly worse, especially when using simpler simulation such as RANS. However, the afore-mentioned agreement is severely reduced when dealing with more complex phenomena such as gust winds. In this case, the data exhibited a significantly worse agreement with the predicted data obtained with the RANS approach. In fact, the authors concluded that RANS was unsuitable for the prediction of gust winds. The overall conclusion of the authors is in high agreement with the one presented by Fernando et al. [60], showing that the complexity of urban flows is the main limiting factor when dealing with urban flows, making it very difficult not to rely on a multi-method approach including in-field full-scales measurements.

Reduced-Scale Measurements
A natural variation of the full-scale experimental approach consists of testing a reduced model in an open environment. On the one hand, we know that urban turbulent flows exhibit a significant degree of complexity and thus are difficult to reproduce in a controlled environment. On the other hand, full-scale measurements tend to involve significant costs and their applicability is typically constrained to the study of the area where they are carried out. Reduced-scale open-environment experiments are an alternative way to approach the experimental study of urban flows. As far as the data-acquisition system is concerned, there is no significant difference, in principle, between the full and reducedscale approaches. The data-acquisition strategy is specifically designed for the experiment in question.
Hirose et al. [66] presented a project to study wind-induced natural ventilation in cities and how it is affected by the surrounding flow. The experimental setup consisted of a 512 cubical blocks matrix where each block had a height of 1.5 m. The total dimensions of the experimental site were 100 × 50 m 2 , and it was oriented such that the wind direction was most often aligned with the length of the site. The data-acquisition system was composed of 700 sonic anemometers equipped with a TR90-T probe. Those were installed at halfheight distance from the top of the block. In addition, acrylic plates with pressure taps were also installed in the northwestern and southeastern faces of the blocks. Using the aforementioned setup two datasets were acquired. From the analytical perspective, the study rests on two major axes. On the one hand, the study of approaching flow conditions, characterised through the statistical descriptions of turbulent flows. On the other hand, the study revolved around the relation between pressure difference and velocity. The authors examined the probability distribution of the velocity at the horizontal wind direction by means of the streamwise velocity magnitude, the standard deviation as well as the velocity range in both southeastern and northwestern winds. As far as the relation between pressure and velocity is concerned, it was examined by means of the pressure coefficient, obtained for every position using the least-squares method.

Closed-Environment Experiments
Closed-environment experiments are, historically, one of the most widely used experimental approaches in fluid mechanics. As it was introduced in Section 3.1, openenvironment experiments gather the complexity of urban flows. However, it is precisely that complexity that limits the applicability range of such methods and thus the generality of the conclusions that can be extracted. It is at this point where controlled-environment experiments in wind tunnels can be helpful. In fact, by controlling the experiment, one can actually study the underlying properties of turbulence and its relation with exogenous parameters.
Wind-tunnel experiments rely on reproducing a given set of flow conditions in a confined space, i.e., a tunnel, in order to study a particular problem. In the case of urban turbulent flows, wind-tunnel experiments are used to study the fundamental mechanisms of turbulence as well as particular applications. Typically, one can distinguish experimental studies by their applicability range, i.e., whether they focus on the analysis of a general phenomenon or on the study of a particular application. In the case of urban flows, this distinction is visible in terms of flow behaviour and applications.
These measurements are usually employed as a verification method. It is very common to use several approaches, e.g. numerical, experimental, modelling, etc. to assess a given problem. Wind-tunnel testing involves inaccuracies, inherent to the measurement methods and the confinement due to the size of the test section. Thus, it is very important to implement the necessary probe corrections [67,68] for the measured quantities, as well as to account for aspects such as three-dimensional effects in the wind-tunnel test section [69]. Vita et al. [65] aim at reproducing the street-level test in full-scale open-environment with a scaled model such that it can be tested in the wind tunnel. Three different types of probes were considered: first, Irwin probes were used to determine the mean wind speed; although these are not the most accurate sensors for this type of measurement, their accuracy was considered to be sufficient their study. The main advantage lies on their omnidirectionality, which eases the installation process since no realignment is needed. Secondly, multi-hole probes such as Cobra probes, were used to measure the incoming wind speed on the top of the model. Multi-hole probes are typically used in high-resolution measurements of turbulent flows. However, those probes are limited by their insensitivity in flows slower than two meters per second as well as their directionality. The third measurement system was hot-wire anemometry, which overcomes the limitations of the two aforementioned systems. Nevertheless, hot-wire anemometers exhibit limitations in terms of spatial resolution and sensitivity to wind direction [70]. To summarize, Irwin probes and hot-wire anemometers were used independently to measure the flow at pedestrian level, i.e., within the model. In addition, a Cobra probe was used to obtain data of the incoming wind speed. Similarly, the work by Gadilhe et al. [64] provides a clear assessment on the verification of measurements systems and models in urban turbulent flows. The work included a comparison of numerical predictions of wind flow with the data obtained from a wind-tunnel experiment. From the experimental point of view, a predictive model was developed to be able to compare both databases. Then a scaled model of the full-scale site was produced to be tested in the wind tunnel under suburban wind conditions.
Wind-tunnel experiments are also employed to study particular underlying effects of urban flows. For instance, Weerasuriya et al. [71] presented a study on the effect of twisted winds at a pedestrian level using a scaled model of the Tsuen Wan street in Hong Kong inside a boundary-layer wind-tunnel. The setup consisted of a series of wooden vanes twisted with a turning table such that the straight-streamlined flow was curved to obtain twisted wind conditions. The twist was set to obtain four distinct cases, i.e., 15 • and 30 • turning clockwise and counter-clockwise. The measurement system included a five-point probing system to obtain data of mean flow, turbulence intensities and yaw angles. At each point, the aforementioned quantities were measured at 12 discrete heights from 10 mm to 1000 mm for a sampling period of 65 s [71].
Other studies use this kind of experiment to particular applications. In the present work, we will focus on analysing two of the most recurrent applications, i.e., ABL and air-quality experiments. The study of ABLs using a wind tunnel typically requires modification of the pre-existing devices [28]. The main objective is to reproduce the conditions representative of an ABL such that wind-engineering applications can be further developed. Hlevca and Degeratu [28] designed a closed-circuit wind tunnel with a 7 m × 5 m × 3 m test section including slotted walls to reduce the blockage effect [28]. When dealing with turbulent-boundary-layer experiments, the distance between the inflow and the test section is critical. It is important to ensure that the boundary layer can properly develop before reaching the test section [72,73]. The authors face precisely that problem as the distance between the inflow and the first measurement point is only 3 m, which was too short for the boundary layer to reach develop sufficiently [28]. Note that increasing the aforementioned distance is typically problematic as the dimensions of the wind tunnel tend to be constrained. Nevertheless, the authors addressed this limitation by installing a device consisting of a grid of horizontal bars with different squared cross sections installed at the entrance of the test chamber [28].
Air-quality studies are a straightforward application of urban fluid mechanics. This area of study can be approached following different paths. On the one hand, one can focus on pollutants and more precisely on their dispersion. On the other hand, some studies approach this problem with a more general perspective, linking fundamental flow behaviour with air quality. In addition, there is a complete field of research dedicated to the study of pollutants from a chemical perspective. The present review focuses on the first two lines of work, leaving the chemical behaviour of pollutants for future works.
The pollutant concentration present in a given environment determines the air quality. From the experimental perspective, this type of study requires the modification of standard wind tunnels such that the pollutant can be properly introduced. Mo and Liu [74] used a water vapor source to introduce particles in the flow. Water vapor was generated with an atomiser to be released as a tracer from the line source [74]. The amount of released vapor-moisture mixture is controlled by the atomizer which is adjusted to match the needs of the different experiments. Note that the pollutant source needs to be properly designed such that the main flow is not significantly disturbed. Mo and Liu [74] actually estimated that the velocity of the flow induced by the line source is around 3%-6% of the freestream velocity. By increasing the distance between the line source and the test section it would be possible to diminish the induced effects. However, increasing the distance between the source and the test area could modify the state of the particles reaching the test section. As far as the measurement system is concerned, it typically consists of the standard equipment with some additional sensors to measure the pollutant concentration. Mo and Liu [74] used a combination of hot-wire anemometry and humidity sensors to measure both flow velocity fields and a conserved scalar, i.e., vapour concentration. Vapour concentration is determined using a combination of two humidity and temperature sensors. Both temperature and humidity are measured at the test region and upstream of the water-vapor source in order to measure reference values. Using the gathered data, relative humidity is computed to be converted into vapor concentration, i.e., specific humidity. An alternative approach was considered by Talluru et al. [75], who released a tracer gas to simulate the contaminant, and measured the pollutant concentration through a photo-ionization detector (PID). In conclusion, one can see that considering pollutant dispersion within wind-tunnel experiment significantly increases the complexity of both the experimental process and the obtained data. In addition, extensive verification has to applied to ensure that the observed phenomena are indeed representative of the physics under study.
Another approach to the study of air quality in urban environments consists of analyzing the fundamental behaviour of the flow. Instead of focusing on the actual drivers of air quality, e.g., pollutant dispersion, the idea is to focus on flow behaviour, which is then characterized using a model. Hertwig et al. [76] employed the aforementioned approach to study the relation between meteorological variability and air quality in urban environments. They combined open-environment observations with wind-tunnel experiments to compare with the prediction of several models [76]. A 20 m × 3.5 m × 1.5 m open-return boundary-layer wind-tunnel was used to test a scaled model of a specific area in south London. The full-scale domain covers 700 m of mixed residential and commercial neighbourhood [76]. Measurements were taken using a two-component laser-Doppler anemometry system. The aforementioned approach provides relevant results regarding the relation between flow behaviour and air quality. However, air quality is introduced with modelling instead of testing in the wind tunnel.

Closed-Environment Experiments with Particle Image Velocimetry
Particle image velocimetry (PIV) is an experimental method to study a certain region within a turbulent flow. Essentially, reflective particles, i.e., seeding, are introduced in the flow and their trajectories are traced. The position of the particles is measured using a camera system which detects them. Note that the seeding has to be properly adjusted, i.e., a sufficient amount of particles have to be introduced without having a significant impact on the flow behaviour. On the other hand, a scaled model of the geometry needs to be produced. Once again, the model needs to be properly set to mimic the flow behaviour around the obstacles under study. The development of the model is typically simple in the case of urban environments since the obstacle array consists mainly of cubic shapes. However, the size of the urban canopies is usually a limiting factor for such experiments. Nevertheless, the PIV technique, when properly implemented, is a very powerful tool to study complex flows.
Monnier et al. [77] presented a research-oriented study that aimed to test the applicability of PIV methodologies in the investigation of flows over an urban-like obstacle array. The authors relied on work of Nagib et al. [78] to design the experimental strategy. The experiment mainly consists of PIV measurements varying the incidence angle of the incoming freestream flow in a scaled urban environment. The measurement system introduced in the wind tunnel consists of an array of three hot wires mounted over a vertical traverse system enabling measurement of the velocity profiles starting from a position close to the lower wall and extending approximately up to one fourth of the domain height wise dimension [77]. As far as seeding is concerned, an atomiser is used to control the size of the particles. Seeding is injected by three inlets at the floor of the facility. Using PIV methodology, the authors verified the correlation between obstacle spacing and flow regimes presented by Oke [18]. Further analysis also brought significant results on the energy exchange between the different boundary-layer regions. In conclusion, the study presented by Monnier et al. [77] shows how the PIV strategy can produce relevant results that deepen the understanding of the fundamental drivers of urban flows.
Ferreira and Ganapathisubramani [79] examined the application of the PIV methodology to the estimation of the pressure distribution in an urban canopy. The authors compared the pressure field derived from in-plane velocity PIV measurements with numerical-simulation data in order to assess the application of PIV measurements to pressure estimations. The geometrical model consisted of two arrays of obstacles, one with uniform height and another one with variable height. The measurements in each of the arrays are obtained using a planar PIV in the streamwise-wall-normal plane [79]. Lasers were installed over the test section to direct the beam towards a set of optics that allowed to redirect the light sheets to predetermined locations. Cameras were mounted on two-dimensional linear translation stages. One of the cameras was tilted to focus on particular planes [79]. The whole system is mounted such that the canopy is fully illuminated. Figure 5 shows a schematic representation of the aforementioned PIV setup.  [79]. Reprinted from Ref. [79], with permission of the publisher (Springer Nature).
As far as the estimation of pressure fields is concerned, the authors use the velocity fields obtained from PIV measurements to run two numerical simulations. Although pressure fields tend to be neglected in several of the studies, from the flow perspective it is important to characterise them properly.The presence of a pressure field in the urban environment modifies the structure of the turbulence, which affects the transport of momentum and thus a proper characterization is very important to the study of the flow. Furthermore, a deeper knowledge of the pressure field helps to assess separation and reattachment locations. Ferreira and Ganapathisubramani [79] also assess the uncertainty of the methodology, with focus on the velocity fields. The main source of error comes from the use of a planar PIV. We know that urban flows are complex and three-dimensional. In this way, the use of a planar measurement system will inevitably introduce errors. The authors used numerical data to evaluate the magnitude of the error, concluding that the uncertainty remained constrained. Furthermore, Ferreira and Ganapathisubramani [79] analysed the error of the numerical methods employed to obtain the pressure estimation. In conclusion, the methodology, despite the uncertainty, is adequately validated [79]. The study in question shows how the PIV methodology can actually be combined with numerical simulations to obtain specific magnitudes of the flow i.e., mean velocity fields, velocity fluctuations, boundary-layer non-dimensional coefficients etc.

Experiments Combining Various Techniques
The work by Vita et al. [65] was aimed at comparing open-environment and windtunnel experiments within the context of urban flows. The authors measured meanflow velocity as well as the corresponding standard deviation. The mean wind-speed data exhibited a good qualitative agreement with the wind-tunnel data, both having the same trend over the measurement positions. On the other hand, the wind-tunnel data was within the standard-deviation range found in the full-scale measurements at several measurement locations. Vita et al. [65] also compared the results with numerical simulations and concluded that the wind-tunnel exhibited deviations with respect to the full-scale results particularly in the recirculation region. Through analysis of their numerical data, Vita et al. [65] found that the velocity streamlines had a rather complex behaviour. In conclusion, the work presented by Vita et al. [65] provides a clear picture on the suitability of the methods available, i.e., numerical methods and wind-tunnel experiments, to study urban flows. The main conclusion presented by the authors suggests the integration of a multi-method approach to deepen the understanding of urban flows.

Numerical Studies
In Figure 6 we show the instantaneous flow around a wall-mounted obstacle [41] for two cases, i.e., with a laminar and a turbulent inflow, which were obtained with a highly resolved simulation method called direct numerical simulation (DNS) using Nek5000 in a spectral-element mesh with a GLL quadrature. In this high-fidelity simulation approach, the governing equations of the flow (i.e., the Navier-Stokes equations) are discretized and solved numerically in a mesh fine enough to capture the smallest turbulent structures. The multi-scale character of the turbulent flow present in urban environments can be observed in Figure 6, where the vortical structures within the wake exhibit a significant range of sizes and energetic content. This high complexity requires the use of appropriate meshes and numerical methods in order to avoid distorting the underlying physics with numerical artifacts. When properly designed high-fidelity simulations are employed, it is possible to obtain a representation of the flow field with the highest level of detail, well beyond what can be achieved with currently available industrial simulations based on Reynoldsaveraged Navier-Stokes (RANS) models (where no turbulent scales are simulated) or coarsely resolved large-eddy simulations (LES, where only the largest turbulent scales are resolved). Although the flow exhibits certain similarities in Figures 6 (left) and (right) due to the fact that the sharp cylinder edges prescribe the separation location, their different inflow conditions lead to a number of visible particularities in each case. The recirculation region on the leeward side of the obstacle is shorter in the turbulent-inflow case, and interestingly the horseshoe vortex forming around the cylinder is much more prominent in the laminar-inflow configuration. This is due to the fact that, as discussed by Vinuesa et al. [41], the incoming turbulence modulates the horseshoe vortex in Figure 6 (left), thus affecting its formation process. Furthermore, as can also be observed in Figure 6, the wake is wider in the turbulent-inflow case close to the obstacle, whereas the laminar-inflow configuration exhibits a wider wake farther downstream. This is consistent with the assessment made by Vinuesa et al. [41] based on the mean flow. A number of studies have investigated the flow in urban environments through numerical simulations [19,80,81], and in most cases the results are validated with wind-tunnel experiments. These numerical studies have documented relevant effects in urban flows which should be studied in even more detail. It is interesting to differentiate the articles simulating urban flows that simply focus on characterizing the flow in urban environments [82,83] (in some cases even using more realistic representations of the shape and distribution of the buildings following the information contained in real planes of specific cities), from the articles considering more realistic situations, although generally using simpler representations of the urban environment (i.e., arrays of regular obstacles). In the latter case, it is possible to identify (among others) some studies comparing several wind-comfort criteria [84], using stratification for calculating surface cooling rates [85] or studying the concentration levels and dispersion of a pollutant release [86].
The characterization of the flow in realistic urban environments is mainly motivated by studying the influence on small-scale fluctuations of atmospheric flow. One important aspect of flow characterization in realistic urban environments is the effect of small-scale fluctuations, which could be produced (among other factors) by the presence of high-rise buildings; it is important to note that in some cases these can produce gusty winds (Monnier et al. [87]). The origin of such turbulence fluctuations could be produced (among other factors) by the presence of high-rise buildings, which in some cases produce gusty winds. Numerical weather prediction models provide a general understanding of the wind system in urban flows, however these tools cannot properly reproduce wind fluctuations, thus requiring more sophisticated models, such as LES, combined with different strategies. Generally, LES is able to characterize the flow in the urban-roughness sublayer, which is the atmospheric layer from the ground to 2-5 times the average height of the building array. This layer is essential to predict the pollutant dispersion, air quality and weather conditions. Regarding the several techniques used to characterize the flow, it is interesting to mention the following works. Nakayama et al. [88,89] performed LES using the Smagorinsky model [90] for the sub-grid scales to simulate a strong wind event in Tokio using a real meteorological setting performing mesoscale simulations. These authors extended the technique by Kataoka and Mizuno [91] to combine turbulent flows from numerical weather-prediction simulations with computational-fluid-dynamics (CFD) models. Giometto et al. [92] performed LES with two closure models for the subgrid scales (SGS), which are the ones that are not resolved in the simulation and need to be modelled. First, the classic Smagorinsky model [90], which relates the SGS terms to the resolved strain-rate tensor. Since this model is known to be over-dissipative in the near-wall regions, Giometto et al. [92] used a wall-damping function (as described in Ref. [93]), which improves the dissipative properties of the Smagorinsky model. Second, the scale-dependent model with Lagrangian averaging of the coefficient, proposed by Bou-Zeid et al. [94], was also used by Giometto et al. [92]. This model improves the performance of the walldamping function, since it uses information from the smallest resolved scales to calculate the model coefficient during the simulation. These authors combined LES with immersed boundary conditions, representing the buildings by a discrete-forcing immersed-boundary method. Giometto et al. [92] compared the results of the simulation with the datasets gathered in the BUBBLE experiment. The BUBBLE experiment is a multi-institutional effort that aimed to acquire data describing the energetics and dispersion processes in the urban boundary layer. A tower was deployed in the "Sperrstrasse" street canyon in Basel, Switzerland to gather data on the velocity components as well as the virtual acoustic temperature [92]. Data were recovered continuously from 1 December 2001 to 15 July 2002 at six different height levels ranging from 3.6 to 31.7 m height. No filtering was applied to the signal to ensure that the energy conservation property is maintained and that the vertical gradients of the state variables computed from the sampled data are properly represented.
The authors show that this formulation is valid to model the turbulence and dispersion in a real urban roughness sublayer, in particular representing the urban city of Basel (Switzerland). Inagaki et al. [95] studied the validity of the outer-layer scaling for turbulence statistics and coherent structures under neutral stratification combining LES with a Smagorinsky model with the lattice-Boltzman method. The authors used a realistic urban geometry modelling the city of Tokio (Japan). They found energetically predominant structures throughout the entire boundary layer, probably driving the flow dynamics, but this information was insufficient to identify the leading mechanism, i.e., whether it was associated to the inner or the outer layer. Finally, Crylls et al. [96] introduced a new method using a LES approach (with periodic boundary conditions) to model the statistically-steady state of non-neutral urban boundary layers. This method freezes transient simulations into a steady state, showing a good performance in both stable and convective conditions.
Regarding the evolution of pollutant dispersion, the review of Belcher [97] reflects the importance of distinguishing between the characteristics of the concentration distribution in the near and far fields, as shown by Theurer et al. [98]. Close to the source, the dispersion depends on the source position and the shape and distribution of the buildings, while in the far field, a Gaussian plume model is sufficiently accurate to model the dispersion (see more details in Refs. [99,100]). Belcher [97] identified three main dispersion processes in the urban environment: (i) the contribution of the shear layers found at the rooftop in the interaction of the flow above and within the canopy; (ii) the topological dispersion, dividing streamlines when the flow impinges on buildings; and (iii) the presence of secondary sources behind the building wakes. Motivated by the lack of information regarding the effect of wind direction and the three-dimensional structure of the dispersion plume in urban areas, Coceal et al. [101,102] performed DNS to compute the statistics of turbulent flows in arrays of obstacles, and Brandford et al. [103] extended that work, also using DNS, studying the evolution of a point-source release of a passive scalar in a urban environment modelled by an array of regular obstacles (cubes). These authors studied the variation of the mean concentration profiles in three different cases, where the angle of the flow direction and the position and number of the point sources was different (three different pressure-gradient directions). They identified six main processes affecting the plume: (i) skewing, (ii) meandering, (iii) detriment by turbulent dispersion, (iv) entrainment and scalar release producing secondary sources, (v) scalar advection downstream and (vi) lateral and topological dispersion due to turbulent fluctuations. To increase the insight into pollution-dispersion mechanisms, Michioka et al. [80] performed LES using dynamic subgrids-scale (SGS) models for the unresolved scales [104,105] to simulate the flow over an array of obstacles with different aspect ratios, showing that the ratio of the turbulent to the advective pollutant flux was almost constant in all the cases, although the turbulent field was the main contributor to the street-canyon ventilation.
The aforementioned work [80,[101][102][103] assumed that turbulent boundary layer was fully developed over large urban areas, with uniform properties. However, more realistic numerical studies take into account that the boundary layer has to adapt to the changing surface roughness, and only a few studies resolve explicitly the roughness transition to simulate the flow [106,107]. Following this research direction, Tomas et al. [108] performed LES in a simplified geometry consisting of over 20 streets. They employed the eddyviscosity model by Vreman [109] to account for the unresolved turbulent stresses. These authors investigate on a rural-to-urban roughness transition the effect of stable stratification on pollutant dispersion, using two lines sources of pollutant. Their results showed the relevance of the vertical momentum fluxes due to turbulence for dispersion of scalars in urban flows. More specifically, they showed that due to the slower internal boundary-layer growth in stable conditions, concentration remains 17% larger. This can be observed in more detail in Figure 7, which shows instantaneous visualizations in the neutral-and stable-stratification cases. Moreover, the advective pollutant flux is still relevant up to a long distance downstream of the rural-to-urban roughness transition. This work was extended by Tomas et al. [110], considering various urban canopies with obstacles using different aspect ratios. These authors studied in depth the influence of such aspect ratio in the pollutant dispersion, surface forces, velocity statistics and boundary-layer depth, trying to identify the dominant mechanisms of pollutant removal from street canyons and their evolution in the transition region. Using a similar eddy-viscosity-based LES approach, Eisma et al. [111] extended the previous work [110], with the aim of assessing the impact of a sound barrier on the urban region. This type of sound barrier is used to reduce the sound levels from highways on neighbouring cities, and the results by Eisma et al. [111] suggest that they may have a significant impact on the pollutant concentration and air quality in the urban area. More specifically, the presence of the plane fence produces a deeper internal boundary later. It is also shown that the mean concentration fields scale with the internal boundary-layer bulk velocity and depth. These results are different from the work by Vinçont et al. [112], who studied the evolution of the passive scalar in the wake of a two-dimensional obstacle, finding the highest concentration in the wake when the line source was located at the wall downstream the obstacle, in contrast to the previous work [111], where the sound barrier was located downstream of the source. Motivated by the need of performing accurate simulations with a reduced computational cost, several authors have developed urban-dispersion models, which are suitable for air-quality and emergency-response applications, although the continuous search of novel approaches that are able to model complex turbulent flow processes is still a research topic of high interest. For instance, Soulhac [113] developed the governing equations which formed the origin of a new operational dispersion model, SIRANE [114]. Using an array of cubes, Hamlyn et al. [115] developed an accurate and very simple model for dispersion, in good agreement with experimental results. Belcher [97] and Belcher et al. [116] showed the explicit dependence of the concentration on the flow parameters and the geometry by the development of a simple analytical model of passive scalar for a regular street network. These authors found that the near-field regime was leading the flow, were the vertical detainment out of the street network was influencing the flow above. Finally, Goulart et al. [117] created a model to predict mean concentration which was validated with DNS results. The DNS modeled the dispersion of a passive scalar from a continuous source release. The authors identified three regions as function of the movement of the vertical fluxes: the near-field, characterized by vertical detainment; the transition region, characterized by vertical detainment and re-entrainment; and the far field, characterized by similar quantities of vertical detainment and re-entrainment. Based on the DNS results, the model uses a formulation based on volume-averaging the governing equations and includes the flow characterization by horizontal advection and vertical detainment and re-entrainment.
Another important aspect of urban-flow simulations is the Reynolds number: a nondimensional parameter representing the ratio of inertial and viscous effects in the flow, which is proportional to the flow speed and the size of the problem, and inversely proportional to the fluid kinematic viscosity. Typical Reynolds numbers in atmospheric boundary layers (ABLs) are several orders of magnitude larger than those of most urban-flow studies based on turbulence models or wind-tunnel experiments [118]. The computational cost of LES and DNS roughly scales with the square and the cube of the Reynolds number respectively [119], which implies that these simulation approaches are prohibitively expensive for full-scale urban-flow simulations. The mean flow in the so-called inertial sublayer has been widely studied in the ABL community, and it is known that (in neutral conditions) it can be described in terms of a logarithmic law [16]. This law includes two parameters representing the surface-roughness and the surface-displacement lengths, which are both characteristic of the urban canopy under study. The logarithmic description of the inertial sublayer led to the idea of using the so-called wall models, which in principle would allow to replace the whole region very close to the wall by a model defining a surface shear stress matching such logarithmic behaviour. This is the cornerstone of most atmospheric models, which avoid resolving the highly computationally expensive scales close to the wall. This wall-model approach, combined with an eddy-viscositybased LES, was used by Cheng and Porté-Agel [107] to study a simplified urban model at high-Reynolds-number conditions. They defined the shear-stress boundary condition based on the Monin-Obukhov similarity theory [120], and their results are very relevant for the improvement of currently existing urban-canopy prediction models. A similar wall-modeled LES was adopted by Giometto et al. [92] to study a real urban geometry in the city of Basel (Switzerland) at high Reynolds number, for which field measurements are available [121]. In their study they adopted a slightly different approach, considering the LES model by Bou-Zeid et al. [94] and the Moeng model [122] for the wall boundary condition. The analysis by Giometto et al. [92] identified turbulent kinetic energy (TKE) production mechanisms above the average building height, as well as the corresponding transport phenomena towards the lower street levels.

Flow Structures and Data-Driven Methods
Coherent flow structures play a critical role in urban environments, since they are responsible for the mechanisms leading to the wind patterns and gusts, thermal distributions and pollutant-concentration fields. Vortical structures are very important structures in the mechanisms of turbulence. Their identification is essential to characterize the physics of the flow. In particular, they are very relevant in urban environments due to the complexity of these flows. Although these structures have been identified and documented in experimental studies, numerical simulations enable sufficient spatial and temporal resolution to evaluate their effect in the dynamic processes taking place in an urban environment. As discussed in Section 4, Vinuesa et al. [41] compared the coherent structures resulting in the flow around a wall-mounted obstacle when the incoming boundary layer was laminar or turbulent, using DNS. As observed in Figure 6, the sharp edges of the obstacle fixes the separation location regardless of the inflow conditions. However, the resulting structures and the topology of the wake significantly differ in both cases. In particular, Vinuesa et al. [41] identified a complex modulation mechanism [125] in the turbulent-inflow case from the incoming boundary layer on the horseshoe vortex forming around the obstacle. This highlights the importance of having well-defined inflow conditions, first for flow reproducibility, and second to ensure an adequate characterization of these types of complex phenomena. It is also important to highlight the need of employing an adequate resolution in the simulations, which may allow to completely characterize elusive phenomena, such as the interaction of the incoming atmospheric boundary layer with an ice sheet [126], a physical process with very important implications in the interpretation of climate models, since it affects the large-scale atmospheric circulation on continental scale [127].
Returning to urban flows, Hunt et al. [128] conducted one of the first studies aimed at identifying the structures around a wall-mounted obstacle, using flow visualization. In Figure 8 we show the schematic representation proposed by Hunt et al. [128], where a total of 4 vortical structures are identified: the horseshoe vortex (A) which forms around the obstacle; the roof vortex (B) and the vortices on the obstacle sides (C), which have a close interaction with the wake; and finally the well-known arch vortex (D) on the leeward side of the obstacle, which will be discussed in detail below. Note that there are a number of methods used to identify coherent structures, and although they do not exhibit significant differences when applied to turbulent flows [129], they may emphasize slightly different phenomena when applied to urban environments [77]. Furthermore, we argue that the use of adequate thresholds may help to discern intense event inherent to turbulent fluctuations and large-scale effects induced by the geometry. A more recent study by Sousa [130] provides a detailed picture of the various structures around a wall-mounted cube using PIV. Sousa [130] used (among others) the so-called swirling strength [131] to obtain timeaveraged structures. Note that the swirling strength λ 2 ci is essentially the square of the imaginary part of the eigenvalues from the velocity-gradient tensor. Further insight into the characteristics of these vortical structures was provided by Becker et al. [132], who considered a series of wall-mounted obstacles with aspect ratios b/h varying from 0.28 to 4. They employed LDA and oil-film visualizations in wind-tunnel experiments where particular care was taken to ensure inflow conditions representative of urban environments. In Figure 9 we show their oil-film visualizations for an obstacle with b/h = 2 and w b /h = 0.29; in this figure, the AOI is progressively increased from 15 • to 60 • . These results indicate the evolution of the arch vortex, which was extensively characterized by Martinuzzi [133] for AOI = 0 • , for increasing incidence angle. The first observation is the dislocation of one of the legs of the vortex, which progressively rotates until reaching the side of the obstacle for AOI = 45 • . Interestingly, an additional increase of AOI to 60 • results in a displacement of the vortex leg to the top of the obstacle, a fact that significantly affects the momentum-transfer mechanisms within an urban canyon. Note that other studies assessing the characteristics of the arch vortex include Refs. [134,135]. As discussed by Monnier et al. [77], there is a wide range of criteria to identify vortical structures, and perhaps the most straighforward approach is to evaluate the vorticity components. In their study, they considered a wind-tunnel test emulating the geometry of the Mock Urban Setting Test (MUST) experiment, which was conducted at the Dugway Proving Ground in Utah as documented by Biltoft [61]. In this simplified urban environment, Monnier et al. [77] considered a set of buildings with b/h = 5 and w b /h = 1, arranged into four rows and three columns. Regarding the building arrangement, they used s/h = 2 and they considered two values of /h, i.e., 1.5 and 4. By visualizing the wall-normal and spanwise mean vorticities, Ω z and Ω y respectively, they could also identify the arch vortex and assess the effect of AOI. Perhaps the most relevant contribution of this visualization is the fact that it highlights the relevance of using a threshold to identify coherent structures. In this case, they considered the modulus of the spatiallyaveraged vorticity vector Ω = Ω 2 x + Ω 2 y + Ω 2 z as a threshold, which allows to clearly see the rotation of the arch vortex when changing AOI. Note that using local thresholds (as opposed to global) can be quite important in order to adequately identify coherent structures, so as to take into account the non-homogeneity of the flow induced by the walls (as it is common practice in wall-bounded turbulence [136][137][138][139]).
Another popular method for vortex identification was proposed by Hunt et al. [140] and it is commonly known as the Q criterion. This method is based on the second invariant of the velocity-gradient tensor Q = −1/2(S 2 + R 2 ), where S and R are the symmetric and antisymmetric parts of the tensor as discussed above. The Q criterion was used by Krajnovic and Davidson [141] to identify the vortical structures around a wall-mounted cube; note that this technique has the advantage of correctly identifying vortices even under high shear. Another popular method for vortex identification is the so-called λ 2 criterion by Jeong and Hussain [22], which is based on analyzing the eigenvalues of S 2 + R 2 . In particular, they define a vortex core as the connected area where S 2 + R 2 has two negative eigenvalues. If the eigenvalues are defined as λ 1 ≥ λ 2 ≥ λ 3 , the vortex core can be identified as the region where λ 2 < 0. There are other vortex-identification methods based on the velocity gradient, such as the ∆ criterion by Chong et al. [142] or the approach by Kida and Miura [143] based on regions of low pressure combined with a swirl condition.
Besides the previously mentioned methods, which are based on the velocity gradient, there is another approach which is based on an integral quantity instead: the normalized angular momentum Γ 1 proposed by Graftieaux et al. [144], which is defined as follows: where D is a small integration domain centered around point x p , u(x) is the velocity vector, × denotes vector product and | · | is the vector modulus. This method was used by Monnier et al. [87] to identify the various relevant regions in the flow between two obstacles, as can be observed in Figure 10. On the leeward side of the upstream obstacle the flow exhibits the arch vortex, with the two legs having flow rotating in the wall-normal direction and the roof exhibiting flow rotating in the spanwise direction. Note that the modulus of the normalized angular momentum Γ 1 is used to define the visualized isosurfaces. The windward side of the downstream obstacle exhibits a region of high spanwise turbulence intensity, and this system is flanked by two regions of high streamwise fluctuations on both sides of the obstacles. Although this flow is observed for an AOI of the incoming ABL of 0 • , see Figure 10 (left), the various elements become reoriented for progressively higher AOI, as observed in Figure 10 (right). This reorientation was also documented in a simpler setting by Rempfer et al. [145]. Note that this particular urban geometry was designed to emulate the flow behaviour in a downtown area, with relatively high buildings compared to the street width, and therefore it exhibits skimming-flow regime [19]. One approach to obtain deeper insight into the flow is by analysing instantaneous three-dimensional fields with data-driven methods, and identifying the coherent structures responsible for the observed behaviour [136][137][138][139]. Furthermore, we believe that a promising approach to gain understanding on the dynamic processes taking place in urban flows is to perform modal analysis, as discussed in the Appendix A. For instance, recent work by Amor et al. [146] has revealed the potential of using dynamic mode decomposition (DMD) to understand the complex physics of the wake in a wall-mounted square cylinder. On the other hand, machine-learning methods have started to receive significant attention in the context of predicting and modelling complex flows. For instance, convolutional neural networks are being used for non-intrusive sensing in turbulent flows, as discussed in detail by Guastoni et al. [147,148]. Convolutional networks are able to exploit the spatial features in the data, in this case the turbulent coherent structures, to perform better predictions. Due to the fact that neural networks are able to perform non-linear predictions, it is possible to use information measured at the wall to sense the flow far from it. The methods based on convolutional networks significantly outperform traditional linear methods, such as linear stochastic estimation (LSE) [149] or linear transfer functions [150], because of the inherent non-linearity of turbulence. Other machine-learning approaches are being developed to perform the modal decomposition of complex flows. For instance, Murata et al. [151] and Eivazi et al. [152] have shown the feasibility of using autoencoders to decompose complex flows into very compact mode representations. Autoencoders can perform nonlinear mode decompositions, as opposed to classical methods, such as proper orthogonal decomposition (POD), where such a decomposition is linear. Although some challenges still exist in terms of the hierarchical relationship of the obtained modes, we believe that autoencoders will really bring in additional insight in terms of the modal decomposition of urban flows. Another aspect required for the accurate prediction and modelling of urban flows is the treatment of time. Srinivasan et al. [153] have shown the possibility of using deep recurrent neural networks, in particular the long short-term (LSTM) network, to predict the temporal dynamics of turbulent flows. Interestingly, this type of network is able to exploit the temporal information in the data, providing excellent predictions of statistical quantities, as well as of the dynamics, in terms of Poincaré maps and Lyapunov exponents. Another interesting approach to model temporal dynamics is based on the Koopman operator with non-linear forcing, as discussed by Eivazi et al. [152]. This approach can provide better short-term predictions than recurrent networks, at a lower training cost. To conclude, we believe that the data-driven approaches discussed above, in terms of non-intrusive sensing, modal decomposition and temporal modelling, may pave the way to more reliable flow predictions in urban environments.

Conclusions and Outlook
In this review we have thoroughly discussed a wide range of available experimental and numerical databases of urban flows, as well as some of the most prominent physical results obtained from them. After discussing all these results, the main conclusion behind this review is that there is still extensive work to do in order to properly develop the research line of urban flows, with the aim to construct sustainable environments in the future. The high complexity of this type of flows and the wide range of possible studies that can be carried out, reflects that this research community has conducted relevant investigations in various research lines, working in parallel, which still need to converge to a single line of action comprising all the findings. This convergence process would unveil the main mechanisms responsible for the dynamics of air pollution, and would enable new strategies to attenuate this undesirable and harmful effect. Among the different open research lines to continue developing, finding accurate predictive models in turbulent flows to properly reproduce pollutant dispersion and developing more robust remotesensing methods would help to achieve new conclusions and new databases that should be analysed in detail. These would provide a more detailed physical insight into the physics responsible for pollutant concentration and thermal distributions. If the flow physics is understood, then it is possible to control it, preventing in this way the harm associated to air pollution.
However, the very large number of spatio-temporal scales defining turbulent flows motivate researchers to search for new tools to identify the main patterns driving the flow dynamics. There are several methods to extract the main features describing the flow [154], although during the recent years the community has paid special attention to novel data-driven techniques to perform modal decompositions of the flow. Combining the acquired physical knowledge of these techniques with new machine-learning strategies it is possible to create fast and efficient reduced-order models (ROMs) accurately modelling the flow [155]. Singular value decomposition (SVD) [156], proper orthogonal decomposition (POD) [13] and dynamic mode decomposition (DMD) [14] are among the most popular techniques generally used for the analysis of various types flows. Their data-driven nature enables their application to a wide range of problems [157][158][159], and they are powerful tools for the identification of flow structures and the construction of ROMs. Their great success has encouraged researchers to develop more robust variants of these methods, exhibiting accurate results in the analysis of complex flows (turbulent, multi-scale, noise). Among the most popular due to their high efficiency in the analysis of turbulent flows it is important to mention spectral proper orthogonal decomposition (SPOD) [160] and higherorder dynamic mode decomposition (HODMD) [161]. These two techniques combine the benefits of POD with the spectral properties of DMD, providing accurate descriptions of the leading modes in turbulent flows driving the deterministic flow dynamics, free of noise and small flow scales [162,163]. More details about the algorithm and the performance of the aforementioned methods is presented in the Appendix A, where also some other variants of POD and DMD are briefly introduced. Based on the understanding of two basic concepts, the spectral and spatial flow complexities, defined in Ref. [161] and adapted to the identification of complex flows in Ref. [163], both SPOD and HODMD are suitable tools for the analysis of urban flows. However, it is important to remark that turbulent flows are composed of both deterministic and stochastic motions, both defining the flow physics. The present methods only identify deterministic motion, which is related to the large-scale coherent structures present in the flow, while the small-scale flow structures, generally understood as low-energy structures, are omitted. Nevertheless, understanding the flow dynamics driving the large size, most energetic, coherent structures present in the flow, would shed light on the physical mechanisms in charge of the pollutant dispersion (as it was tested in other complex problems in [164][165][166]). Hence, the new future trends should focus on continuing to apply and develop highly efficient data-driven tools to identify the main flow patterns in urban flows, not only focusing on the aforementioned classical tools, but also exploring new algorithms of machine learning-for instance, those based on nonlinear functions. Combining the fields of artificial intelligence and fluid dynamics could provide new models, with reduced degrees of freedom and with predictive capacities, that would also consider stochastic flow motions, providing additional insight into the main mechanisms related to urban pollution, but this field is still barely unexplored. The great success of artificial intelligence in several other different fields (i.e., banking, earthquake predictions and computer vision, to name a few), suggests the high potential of this new research line which should be investigated in the near future. Identifying the main flow patterns with high accuracy, and developing high-fidelity ROMs, will shed light into new mechanisms related to the flow behaviour, thermal effects and pollutant dispersion in this type of environment.  and efficient reduced-order models (ROMs). Using data-driven techniques to extract the main features driving the flow, and combining the acquired physical knowledge with new machine-learning strategies to create ROMs is a research topic of high interest that should be explored more in detail. This Appendix presents a brief summary of some relevant datadriven techniques based on modal decompositions exhibiting good performance when analyzing complex flows [157,158,160,162,163,167] (see more details and other techniques in Refs. [155,159,168]). These methods are powerful tools that could be potentially used in the construction of new ROMs for urban environments, among others applications.
For convenience, to apply the following data-driven tools, the data are organized into a snapshot matrix, which considers K snapshots evolving in time as: where v k corresponds to the field variable evaluated at time instant t k , for simplicity defined as v k = v(t k ). The dimension of the snapshot matrix is J × K, where J represents the spatial degrees of freedom of the analysed data and K is the number of snapshots, as defined before. Singular value decomposition (SVD) [156] and proper orthogonal decomposition (POD) [13] are well-known techniques that are generally used to create reduced-order models [169] and to identify flow structures and patterns [170,171]. The main characteristic of these two techniques are their ability to represent high-dimensional solutions through low-dimensional models, defined by a few SVD/POD modes, which are associated with a high energy content as detailed below.
POD decomposes the flow field v(x, y, z, t) (in fluid dynamics generally the velocity vector) into temporal, c j (t), and spatial orthogonal modes, Φ j (x, y, z), as: where the selection of such spatial modes is optimal in terms of the system energy (kinetic energy if the flow field under consideration is the velocity vector). POD modes can be identified using two types of algorithms, the classical POD method [13], which is based on the covariance of a state vector changing in time, and the snapshot method [156], which is similar to the SVD method. Hence, while POD is the method providing the expansion (A2), SVD can be considered as one of the algorithms providing the POD modes.
Both the classical POD algorithm and the SVD method are based on the definition of the covariance matrix as C = X X (with dimension J × J), where (·) denotes the matrix transpose and X is the snapshot matrix defined in (A1). By definition, this matrix is symmetric and positive-definite, hence its corresponding eigenvalues and eigenvectors will be real numbers, and the eigenvectors will be orthogonal among them, forming an orthogonal basis of modes.
The classical POD method solves the following eigenvalue problem: where the optimal orthogonal modes Φ j (POD modes) that better approximate a dynamic field, are related with the largest eigenvalues λ j . In three-dimensional turbulent flows, this method is prohibitively expensive since it is based on solving the eigenvalue problem of the covariance of a state vector changing in time, with dimensions proportional to the spatial degrees of freedom (J × J). Hence, the snapshot method should be used instead. This method is based on the fact that the most energetic POD modes are the same as those obtained with the transposed covariance matrix as : The dimensions of C is proportional to the snapshot number from Equation (A1), i.e., K × K (let us remember that the dimension of X is J × K), in fluid dynamics generally being K << J. Moreover, when the snapshot matrix (A1) is composed by velocity vectors, this covarianze matrix contains the variations of the square of the velocity vector, which agrees with the definition of kinetic energy. In such cases, the POD modes are ranked by levels of kinetic energy.
SVD is a type of factorization that captures the directions of a matrix in which vectors can either grow or shrink. These directions are given by the eigenvalues and eigenvectors of a rectangular matrix. Similar to the snapshot method, SVD applied to the snapshot matrix (A1), decomposes the flow field into spatial (W) and temporal (T) modes and singular values Σ, as: where W W = T T = I, with I being the N × N unit matrix, and Σ is a diagonal matrix composed by the singular values (σ 1 , · · · , σ K ), which are ranked in decreasing order. This method is strongly connected to the previous eigenvalue problems, where σ 2 j correspond to λ j , the columns of W are the POD modes Φ j (ranked consistently with their corresponding eigenvalues) and the columns of T are the temporal modes Ψ j .
The number of SVD (or POD) modes most relevant to describe the flow can be identified using several criteria. Several methodologies are described in Ref. [172], where SVD and POD methods are also known as principal component analysis (PCA). When solving fluid dynamics problems the standard SVD error, estimated for a certain tolerance ε 1 , defines the most relevant modes in POD analyses based on the singular values as: The most relevant modes can then be used to provide low-dimensional approximations of complex flows.
The good performance of POD/SVD algorithms has been shown in the identification of coherent structures in a wide range of applications modelling complex flows [169,170]. However, motivated by improving the accuracy of the method to identify patterns in turbulent flows, some authors have developed several variants of the POD algorithm. To name a few: the multi-resolution POD algorithm [173], where before proceeding with the POD analysis of the data, the correlation matrix C is split into the contribution of different flow scales using a filter (preserving symmetry of the correlation) or a Wavelet Transform; the balanced POD [174,175], which is generally used to produced high-fidelity reduced systems combining POD with the balanced truncation in Control Theory [176]; or the spectral proper orthogonal decomposition [160], which is explained more in detail in the following section.

Appendix A.2. Spectral Proper Orthogonal Decomposition
Spectral proper orthogonal decomposition (SPOD) [160] is an extension of POD recently introduced for analysis and pattern identification in turbulent flows [162,[177][178][179]. This method calculates the POD modes in the frequency domain, combining the good properties of POD with the information contained in a temporal modal expansion as in dynamic mode decomposition (described below).
The algorithm starts from the snapshot matrix (A1), with dimension J × K (J represents the degrees of freedom of the computational domain and K is the number of snapshots), which is divided into several windows with dimension J × K , where K < K. Applying the Fast Fourier Transform (FFT) to each one of such windows, the main frequencies are calculated. A dominant frequency is selected from the multiple solutions in the different windows and the associated modes are used to create a new snapshot matrix. Finally, POD is applied to the new snapshot matrix, reducing the complexity of the data and removing the noise and the spatial redundancies. This method provides low-rank approximations of the solution, free of noise, small flow scales or spurious results. However, application of FFT to reduced groups of data can be computationally expensive, particularly in threedimensional computational domains where J may be large. A good option could be inverting the order of the steps performed in algorithm. In other words, it is possible to reduce the degrees of freedom from J to N by first applying the POD method. This is the essence of the higher-order dynamic mode decomposition (HODMD) algorithm presented in the following section. Dynamic mode decomposition (DMD) [14] is a technique recently introduced to identify flow structures. The method decomposes spatio-temporal data, v(x, y, z, t k ), into Fourier-like modes as: v(x, y, z, t k ) M ∑ m=1 a m u m (x, y, z)e (δ m +iω m )t k , for k = 1, · · · , K, where u m (x, y, z) are the DMD modes weighted with the amplitudes a m , ω m are their corresponding frequencies and δ m are their associated growth rates, representing the temporal growth or decay of the modes in time. The DMD algorithm introduced by Schmid [14] is based on the linear relationship between two consecutive snapshots via the linear Koopman operator R. Starting from the snapshot matrix (A1), where the snapshots are organized equidistant in time with time step ∆t, the DMD method is defined using the following Koopman assumption: where the matrices V K−1 1 and V K 2 contain from the first to the last but one snapshots, and from the second to the last snapshots of the data base, respectively. The Koopman matrix R contains the dynamics of the system, while the DMD frequencies ω m and growth rates δ m are the computed eigenvalues and the computed eigenvectors are used to construct DMD modes u m (x, y, z). Recently, Le Clainche & Vega [161] introduced the higher-order dynamic mode decomposition (HODMD) method, which is an extension of DMD for the analysis of complex flows (turbulence, multi-scale, flows in transitional regime,...) and noisy experimental data. This method combines the Koopman assumption (A8) with the Taknes' delayed embedded theorem [180]. Hence, HODMD relates d time-delayed snapshots (sub-matrices) using the high-order Koopman assumption as: which can be understood as the window-shift process carried out in power-spectral density (PSD). Finally, as in DMD, the dynamics of the system is contained in the several Koopman operators, R 1 , · · · , R d , the eigenvalues of which represent the DMD frequencies and growth rates and their eigenvectors are used to construct the DMD modes, as in standard DMD. The HODMD algorithm proceeds in two steps (see more details in Ref. [161]): • Step 1: dimension reduction. Applying truncated singular value decomposition (SVD) to the snapshot matrix V K 1 yields: The matrixT K 1 is the dimension-reduced snapshot matrix. The number of SVD modes retained in this approximation N is defined as the spatial complexity. These modes are selected as in the SVD algorithm presented in Appendix A.1. A (tunable) tolerance ε 1 estimates the standard SVD error, as described in Equation (A6). • Step 2: the DMD-d algorithm for the dimension-reduced snapshots. The highorder Koopman assumption is applied to the reduced snapshot matrix, resulting in: After some calculations, the several Koopman operators R 1 , · · · , R K are grouped into a single matrix, the eigenvalue problem of which provides the DMD modes, frequencies and growth rates that define the DMD expansion (A7). This expansion is sorted in decreasing order of the mode amplitudes and it is further truncated by eliminating the modes such that: for some tunable parameter ε 2 . The number of modes retained in this expansion, M, is called as the spectral complexity.
When d = 1 in equation (A11), HODMD reduces to standard DMD, defined as: while for d > 1, HODMD can be understood as the result of applying standard DMD to a set of enlarged snapshots, defined by the delayed snapshots. Hence, HODMD combines the advantages of standard DMD with some consequences of the delayed-embedding theorem by Takens [180]. The spatial and spectral complexities determine the performance of the HODMD algorithm. When the spectral complexity is larger than the spatial complexity, M > N, DMD fails (see additional details and examples in Ref. [161]), thus requiring to use the DMD-d algorithm, with d > 1, to overcome the problem. In complex flows (turbulent, noisy, etc.), the parameter ε 1 is selected so as to retain only the large scales in the flow. Hence, at this step the spatial complexity N can be smaller than the spectral complexity M, which is generally very large when modelling complex dynamics (see more details in Ref. [166]). Consequently, HODMD is a suitable tool to identify flow structures in noisy experimental data [164], flows in transitional regime [181] and in other complex flows [166,182,183], even with non-equidistant data in time [163], which is one of the limitations of classical DMD methods.
Finally, it is important to mention the spatio-temporal Koopman decomposition (STKD) [184], an extension of HODMD which decomposes the flow as: v(x, y, z s , t k ) M,B ∑ m,b=1 a mb u mb (x, y)e (δ mb +iω mb )t k +(ν mb +iβ mb )z s , for k = 1, · · · , K, and s = 1, · · · , S, with ν mb and β mb being the spatial growth rates and wavenumbers associated to the spanwise direction. Similarly, the flow can be decomposed along the remaining spatial directions of the analysed data as: v(x i , y l , z s , t k ) M,B,P,Q ∑ m,b,p,q=1 a mbpq u mbpq e (δ mbpq +iω mbpq )t k +(ν mbpq +iβ mbpq )z s +(µ mbpq +iγ mbpq )y l +(σ mbpq +iα mbpq )x i , for k = 1, · · · , K, s = 1, · · · , S, l = 1, · · · , L and i = 1, · · · , I, with α mbpq , γ mbpq and β mbpq being the spatial wavenumbers and σ mbpq , µ mbpq and ν mbpq the corresponding spatial growth rates. Details of the STKD algorithm can be found in Ref. [184]. The essence of the method is to apply the high-order Koopman assumption presented in Equation (A11) to the spatial modes (weighted with the singular values and properly re-organized) of the SVD in (A10) from Step 1 in the HODMD algorithm. A novel Koopman operator defines a new eigenvalue problem, the eigenvalues of which represent the spatial growth rates and wavenumbers in the previous expansion. This method has been successfully used to identify the main flow structures in complex flows [163,166] and to identify the flow instabilities leading the flow transition in the wake of a wind turbine [165].

A.4. Some Variants of DMD Not Related to Time-Lagged Snapshots
In recent years, DMD has been used in a wide range of applications in fluid dynamics to identify coherent structures [157] and global instabilities [185,186]. The high potential of this tool has motivated several researchers to define new algorithms which widen the type of applications of the tool, providing in some cases more accurate solutions. To mention a few of them: • Sparsity-promoting DMD [187]. This method uses convex optimization techniques to identify a smaller set of important modes. • Extended DMD [188]. This algorithm includes more basis functions in the standard DMD approximation, a fact that allows the method to retain more modes, thus enabling the description of more complex dynamical systems. • Optimized DMD [189]. An optimization problem is solved to compute the DMD expansion (A7). Both DMD and HODMD are purely linear-algebra-driven approaches compared to optimized DMD. Furthermore, other authors [190] have developed other linear-algebra-oriented algorithms which treat separable nonlinear least-squares problems. • DMD variants to treat noisy data. Several authors have put an effort on identifying and removing noise from the analysed data. Among other techniques it is important to mention the method by Dawson et al. [191], which characterizes the noise properties; the total least-squares DMD by Hemati et al. [192], which combines standard DMD with total least squares; and the method by Takeishi et al. [193], which combines standard DMD with a Bayesian formulation. • Multi-resolution DMD (mrDMD) [194]. The main idea behind mrDMD is to separate the high-and low-frequency events which generally occur in complex flows. The algorithm divides the snapshot matrix (A1) into several segments to identify highand low-frequency modes depending on the segment length. The resulting DMD expansion (A7) is then represented by the various sub-expansions of DMD modes with different frequency ranges. A similar idea is behind the multi-resolution POD algorithm [173].