### #4 POLLUTION

Elaboration of predictive models for the propagation of liquid and atmospheric pollutants in the city of Lisbon.

**Objective**

To develop predictive models for the propagation of liquid and atmospheric pollutants.

**Research question**

How is the propagation of air and liquid pollutants in the city, in case of an accident with hazard substances?

**Business understanding**

For the definition of the #4 Pollution use cases, were done meetings with Civil Protection Services from the Lisbon City Council (SMPC) and was discussed the importance of pollutants dispersion case of accidents.

For this case two possible situations that SMPC considered relevant were for one side the vulnerability of an hospital and population in Alcântara area (Lisbon) in case of a leakage of natural gas, and the population vulnerability that lives near some segments of the railroad, considering that there’s a substantial train traffic that transports hazardous materials, crossing the city.

Taking into consideration these aspects were defined two use cases: 1) the simulation of the propagation of natural gas in Alcântara area after a disruption in the gas network; and 2) the simulation of the propagation of nitrobenzene in Roma-Areeiro train station after a train accident.

The Barcelona Supercomputing Center has been in charge of the #4 Pollution use cases, being the primary goal of this part of the project the elaboration of predictive models for the propagation of liquid and atmospheric pollutants in different areas of Lisbon.

The ultimate objective of these predictive tools is to simulate different pollutant emission scenarios allowing the creation of emergency measures to minimize the impact of an eventual hazard on public health and other citizen's life aspects.

**Data understanding**

For the development of the #4 Pollution use cases was requested data to the Lisbon City council and to the Portuguese Institute for Sea and Atmosphere (IPMA) (Table 31).

###### Table 31. Datasets necessary for the implementation of the analytical model for #4 Pollution use cases

Use case 1: Simulation of air pollutants in the city of Lisbon after a rupture in natural gas network

This case study is intended to model a leakage of natural gas in Alcantara (with residential, tourist, harbour area and a hospital), with a temporal frequency of 1 minute to 15 minutes, in an area of 1 km2 calculated from the location of the incident (Figure 35).

###### Figure 35. Study area of the use case in which was simulated a leakage of natural gas after a disruption in the gas network (Alcântara)

Some assumptions have been made to determine the necessary boundary conditions for the case study modelling. In this case, the assumptions were:

- A standard natural gas pipe for urban areas was considered;
- Pipe diameter 0.2 m;
- Pipe distribution pressure: 200 MPa.

These assumptions allowed us to determine the mass flow rate of the gas leakage (Lu, Zhang, Yan, Li, & Zhao, 2014):

Q=15000 m3/h

Considering the pipe diameter, this entails an inflow velocity of purely natural gas of:

u=132 m/s

The following coordinates projected in the ETRS89 / Portugal TM06 (EPSG: 3763) coordinate system, were considered as the source of the incident:

Y: -106559,43952

X: -90388,247875

Use case 2: Simulation of liquid pollutants in the city of Lisbon after a train accident with hazardous liquids

Most hazardous materials circulate in Lisbon by train, with their departure and destination in Alcântara station. During their route, the trains pass by the stations of Oriente, Braço de Prata, Roma Areeiro, Entrecampos, Sete Rios, Campolide and Alcântara. Nitrobenzene is one of the transported substances, among others. In this use case is intended to model the propagation of nitrobenzene after a train accident in Roma-Areeiro station with a temporal frequency of 1 minute to 15 minutes, in an area of 1 km2 calculated from the location of the incident (Figure 36).

###### Figure 36. Study area of the use case in which was simulated a leakage of nitrobenzene after a train accident in Roma-Areeiro train station

The initial definition of the case left a broad margin for interpretation.

Therefore, several assumptions have been made to obtain the required input data to perform the numerical modelling, which is the pollutant emission rate to the atmosphere and its geometrical distribution. These are the scenario conditions that have been assumed:

- The crashed train is an ISO tank container of 22m3 of capacity;
- The spilled liquid will be spread along the track area;
- The spread fluid film will be around 2 cm thick.

These hypotheses allow the evaluation of the spill evaporation surface, which is assumed to be a rectangular surface aligned with the tracks with a total area of 1000m2.

On the other hand, to determine de evaporation mass flow rate, the conclusions of an experimental investigation with nitrobenzene have been used (Hine, 1924).

The experiment evaluates different scenarios with varying wind speeds and temperature conditions to determine the evaporation mass flow rate for a circular surface of 0.6m of diameter.

Our study used the worst scenario, which is the fastest pollutant release ratio to the atmosphere that is compatible with realistic weather conditions. These conditions are:

- wind speed 5.18 [m/s];
- Avg. air temperature 24.2 [ºC];
- Avg. liquid temperature 25.5 [ºC];
- mean vapor pressure 48.12 [Pa];
- Relative humidity 48 [%].

This leads to a mass flow rate of M=0.445 [mol/h].

Considering that the relative density of the Nitrobenzene saturated vapor at 20ºC is 6.6 (i.e., 1 is the air density), we can derive the volumetric flow rate of substance released in a 1 m2 of surface:

Q = 6.8897 1e-6 m3/s

This volumetric flow allows prescribing the pollutant inflow boundary condition for the dispersion scenario: inflow velocity (u=6.8897e-1 m/s) and pollutant concentration (c=1).

The following coordinates projected in the ETRS89 / Portugal TM06 (EPSG: 3763) coordinate system, were considered as the source of the incident:

Y: -101946,9944

X: -87149,50364

Lisbon meteorological data

Since the pollutant dispersion patterns are mainly driven by the local micrometeorology, it is essential having insight into Lisbon's predominant meteorological phenomena. To that end, a spreadsheet containing the temporal measurements of three meteorological stations in Lisboa was provided by Lisbon City Council at the beginning of the project.

Specifically, the measurements were retrieved from three IPMA weather stations of IPMA (Ajuda, Geofisico, and Gago Coutinho). The measurements are given in hourly time series from 01/01/2013 to 31/12/2018 for each station. The following data was contained in each measurement in the dataset:

- Measurement time and date;
- Station id;
- Station name;
- Temperature (oC);
- Humidity (%);
- Precipitation (mm);
- Wind speed (m/s);
- Wind direction (o).

This data allowed to determine the dominant wind directions and intensities at each station which were used to model the different use case scenarios. For this purpose, the wind roses have been generated, plotting the statistical frequency of each wind direction and velocity for Ajuda (Figure 37), Geofísico (Figure 38), and Gago Coutinho (Figure 39) weather stations.

###### Figure 37. Ajuda station wind rose., along with wind speed (m/s) for each direction

###### Figure 38. Geofísico station wind rose., along with wind speed (m/s) for each direction

###### Figure 39. Gago Coutinho station wind rose., along with wind speed (m/s) for each direction

All the meteorological stations show a very similar behaviour regarding the wind directions, being the east direction the dominant one. This direction matches Tagus river direction. Thus, all the simulated scenarios will be run with this wind direction since it is the most likely to happen in a particular event.

Data preparation

Having the urban geometry and topography is essential to reproduce the use case scenarios in realistic conditions. To that end, Lisbon City Council, provided CAD data of the two areas of study. The data contained the building geometry and the topography separately and at different scales and reference origins. The data was adequately processed to match both domains (Figure 40).

###### Figure 40. Original CAD data processing for the Alcântara use case

After performing the preliminary simulations to assess the modelling strategy performance, it was concluded that upstream geometry was necessary to allow the boundary layer to develop a realistic urban behaviour before reaching the study area.

For this purpose, extended data was required for both use cases. Particular emphasis was placed on the upstream region. The final geometry for Alcântara and Roma-Areeiro is represented in Figure 41 and Figure 42, respectively.

###### Figure 41. Alcântara area final geometry. The original area of study is represented by the red square

###### Figure 42. Roma-Areeiro area final geometry. The original area of study is represented by the red square

Geometry meshing

The final step before running the simulation is to create the simulation domain and the computational mesh on which the discrete solution of the problem will be obtained.

Using adequate domain size and grid resolution is essential to obtain correct solutions. To that end, we followed the requirements established in the AIJ guidelines for practical applications of Computational Fluid Dynamics (CFD) in urban environments (Tominaga et al., 2008).

According to these guidelines, the following geometrical and mesh parameters have been used:

- Domain size: 10000x10000x1000 m;
- Grid size at the pedestrian level: 1m.

The horizontal domain sizes are ten times the size of the study area. These sizes are significantly larger than those required in the AIJ guidelines (Tominaga et al., 2008). However, we prescribed them to allow the use of periodic conditions in the wall-parallel directions. This strategy will be described in more detail in the modelling section.

In Figure 43 is represented the computational mesh of the Alcântara area and a detail of that mesh in Figure 44.

###### Figure 43. Computational mesh of the Alcântara area

###### Figure 44. Detail of the computational mesh of the Alcântara area

**Modelling**

This part of the project aims to create an innovative, accurate and reliable modelling framework to predict the dissemination of atmospheric pollutants in urban environments. To that end, CFD will be used as the primary modelling strategy to meet the accuracy requirements. CFD has been traditionally used to perform microscale meteorological predictions.

By microscale, we mean spatial domains of the kilometer order with a computational resolution that can capture street-level phenomena.

The numerical meteorological predictions include at least the wind velocity and the pressure fields within the computational domain. The data obtained from these simulations can be of capital interest to understand the flow dynamics in urban and suburban areas since obtaining experimental measures with the same degree of detail and resolution is unfeasible.

On the other hand, microscale meteorological simulations can be used as a transport field to predict pollutant dispersion patterns, among other applications. Nevertheless, within the CFD category, there is a vast range of different methodologies with different degrees of accuracy and reliability that can be used.

In general, the most important selection criteria is the availability of computational resources.

Traditionally, most of the works related to microscale meteorological prediction in the literature use Reynolds Averaged Naver Stokes (RANS) simulations to perform the flow predictions given their relatively low computational cost (Li et al., 2010). However, the RANS methodology introduces strong physical assumptions in its formulation, and thus, this choice comes at the expense of the simulation's accuracy and reliability.

The continuous increase in computational power availability and efficiency allows higher fidelity simulation techniques for microscale meteorology prediction in urban environments nowadays. Specifically, Large Eddy Simulation (LES) techniques produce numerical solutions with high resolution in space and time. In LES, most relevant physics are explicitly resolved, introducing only weak physical assumptions in the model. This makes the solution obtained through LES much more accurate and reliable than those computed with RANS models.

According to Toparlar et al. (2017) only 6 out of the 61 most relevant works published on micrometeorology of urban flows use high-fidelity LES simulations, and the analysed computational domain in most cases is restricted to small urban areas (Li et al., 2010; Yaghoobian et al., 2014). To our best knowledge, only García-Sánchez et al. (2018) and Kurppa et al. (2018) have presented fully real urban areas simulated through LES. Precisely, García-Sánchez et al. (2018) reproduced the Joint Urban 2003 experiment, in which actual measurements of wind fields were taken in Oklahoma City. On the other hand Kurppa et al. (2018) resolved the wind field and non-reactive gas pollutant dispersion in an 8.2 km2 area of Helsinki.

For this specific project, we wanted to go further in using high-fidelity CFD techniques in urban area flow simulations by proposing an innovative approach, which covers from the physical and numerical method point of view to practical implementation issues. Given the vast extension of solid-wall areas of these applications, we proposed using wall-modelled Large Eddy Simulation as a methodology to perform the flow predictions.

LES produces high-resolution solutions both in space and time. Since the real physics are well resolved in all the problem dimensions, only weak physical assumptions are needed to close the model, enhancing the reliability of the obtained solutions. Unfortunately, the required spatial and temporal resolutions may entail an unaffordable computational cost, especially at high Reynolds numbers and with the presence of sizeable solid wall areas as in the present study cases. Boundary layers are complex flow structures that are engendered in the vicinity of the solid wall.

The explicit resolution of these flow regions for city-sized computational domains is unfeasible, and alternative strategies had to be used since their effects on the flow are crucial and cannot be neglected. In that sense, wall modelling for LES (WMLES) is a suitable approach to overcome such a significant difficulty.

The main advantage of this strategy is that we can keep the reliability and accuracy of a LES prediction while circumventing the costs of the explicit wall resolution. In Figure 45 is presented the isometric view of WMLES simulation in the Alcântara area. An horizontal section of WMLES velocity at pedestrian level for the same study area is presented in Figure 46.

###### Figure 45. Isometric view of a WMLES simulation in the area of Alcântara. Contours of the velocity magnitude are plotted

###### Figure 46. Horizontal section of the WMLES velocity field at pedestrian level in Alcântara area

In that sense, for this project, we used Alya, an existing in-house HPC code (Vázquez et al., 2016) developed at CASE Department. The resolution of wall-bounded turbulent flows through WMLES has been extensively investigated by the members of the Alya development team, and several solutions have been proposed, implemented, and validated in the code.

Turbulent flows are the main physics involved in the proposed simulations. The extreme spatial complexity and unsteadiness of turbulence make the governing equations' numerical solution (filtered Navier-Stokes equations) highly sensitive to the discretization method. Spatial discretization schemes may introduce numerical diffusivity to the solution, a completely unphysical phenomenon that severely distorts the solution. Alya features low-dissipation schemes for the finite element method (Lehmkuhl, Houzeaux, Owen, Chrysokentis, & Rodriguez, 2019), in which the code is based. This approach is crucial for the correct modelling of resolved turbulence. On the other hand, Alya also features subgrid (LES) models developed by the members of the BSC team (Lehmkuhl, Piomelli, & Houzeaux, 2019), which, together with the appropriate discretization methods commented above, ensures the code ability in reproducing the flow physics successfully. This makes Alya a reliable code for high-fidelity CFD simulations (Lehmkuhl, Piomelli, et al., 2019; Pastrana, Cajas, Lehmkuhl, Rodríguez, & Houzeaux, 2018; Rodriguez et al., 2019).

Regarding the modelization of the wall region, the main difficulty arises from the extraordinarily complex and tiny size of the near-wall flow structures compared to the whole domain length scale while having an essential influence on the global flow behavior. Wall modeling allows considering their effects in the LES computation through a physical model instead of explicitly resolving them. This introduces a new layer of modelization in the overall physical model, which may undermine the flow prediction's accuracy. The members of the Alya development team have broad experience in developing and implementing wall modeling approaches (Calafell, Trias, Lehmkuhl, & Oliva, 2019) for finite element-based applications (Owen et al., 2020), which minimized the potential effects of this additional source of errors. The models implemented in Alya have proven their reliability when simulating flows similar to those of the current proposal (Calafell et al., 2019; Owen et al., 2020).

Another difficulty in the specific case of time-resolved environmental simulations is the prescription of appropriate time-resolved boundary conditions at the lateral boundaries. In general, the boundary conditions for the microscale domain are provided by external codes such as the Weather Research and Forecasting (WRF) mesoscale model. This model provides time-varying pressure and velocity distributions to be interpolated at the microscale domain boundary (Liu, Miao, Zhang, Cui, & Zhang, 2012). Nonetheless, the characteristic times of the time-dependent variables provided by WRF do not necessarily match the turbulent time scales that should be present in an LES-resolved atmospheric boundary layer (ABL). To overcome these difficulties, we propose an easy-to-implement method based on periodic boundary conditions in all wall-parallel directions. The wind motion is enforced through the geostrophic pressure gradient provided by WRF, which is imposed as a source term to the model. This creates a natural and straightforward coupling between the mesoscale and microscale simulations through a single vectorial value. At the same time, the proper resolution of the turbulent time scales throughout the computational domain is ensured by this method since this approach is widely used for DNS and LES of canonical wall-bounded turbulent flows (Lozano-Durán & Jiménez, 2014).

Physical and mathematical model

The physical and mathematical model implemented in Alya is based on the numerical resolution of the spatially-filtered incompressible Navier-Stokes Equations:

where *u*¯ is the filtered velocity field, *p*¯ is the filtered pressure, ? is the fluid molecular viscosity, and ?sgs stands for the subgrid viscosity.

This is a suitable model to characterize the ABL physics (García-Sánchez, van Beeck, & Gorlé, 2018; Kurppa et al., 2018; Owen et al., 2020).

To close the formulation, a proper subgrid model was used to provide the ?sgs. For the present modelling approach, the Vreman model was used given its proven reliability in modelling atmospheric boundary layers (Owen et al., 2020):

where *c* is a model constant while ?m is the grid characteristic length-scale.

On the other hand, to deal with the difficult near-wall physics, a wall model based on an equilibrium atmospheric boundary layer profile with roughness is used. This approach has proven to be suitable for atmospheric boundary layers (Owen et al., 2020) and urban flows (García-Sánchez et al., 2018).

being *ux *the wall-parallel velocity component, *y *stands for the wall distance, u? is the skin friction velocity, ? is the Von Kármán constant, and *y*0 is the wall roughness height.

Finally, the transport of the pollutant substance within the atmospheric flow has been modeled as a passive scalar through a filtered advection-diffusion equation:

where *c*¯ is the relative substance concentration (c¯=1 represents pure substance), *D *is the diffusivity, *D*sgs is the subgrid diffusivity due to the unresolved scales, and *S *stands for any substance source.

To determine the diffusivity coefficients, a constant Schmidt number has been applied. Again, this approach has been validated to model the transport of pollutant substances in urban environments (Gousseau, Blocken, Stathopoulos, & van Heijst, 2015).

To summarize, the variables that will characterize the micrometeorology and the dispersion problem are the wind velocity field, the pressure, and the pollutant substance concentration at each point of the high-resolution computational grid.

To deal with the complex issue of the boundary conditions for a time-resolved turbulent flow, a novel strategy has been developed for the present simulations. The approach is based on prescribing periodic boundary conditions in all domain’s wall-parallel directions. The wind motion is enforced through the geostrophic pressure gradient provided by WRF, which is imposed as a source term to the model.

This creates a natural and straightforward coupling between the mesoscale and microscale simulations through a single vectorial value (Figure 47). At the same time, the proper resolution of the turbulent time scales throughout the computational domain is ensured by this method since this approach is widely used for DNS and LES of canonical wall-bounded turbulent flows (Lozano-Durán & Jiménez 2014).

###### Figure 47. Wind flow enforcement scheme

Model validation

An a priori validation test has been carried out to ensure that the global modelling strategy performed successfully. The main objective of the test was to ensure that the periodic strategy combined with the physical and numerical model was able to reproduce realistic velocity profiles upstream of the urban area.

To achieve this, the downstream perturbed atmospheric boundary layer has to recover the original logarithmic profile along the periodic loop. For that reason, the urban area of Shinjuku (Tokyo) has been used to validate the strategy (Figure 48).

Its geometry is particularly challenging for the model due to the height of the buildings present in this area, severely disrupting the boundary layer equilibrium profile (Figure 49).

###### Figure 48. Shinjuku (Tokyo) urban area geometry

###### Figure 49. Shinjuku urban area ABL. Top: mildly perturbed ABL featuring small-building urban geometries. Bottom: Strongly distorted ABL due to the presence of tall buildings

Specifically, three domain sizes were tested, 3L, 5L, and 10L being L the urban area horizontal length (1Km). On the other hand, the wall model roughness was set at y0=0.1 m. Below, the mean velocity profiles upstream and downstream of the urban area are displayed for the 5L and 10L tests (Figure 50).

###### Figure 50. left: downstream average velocity profile, middle: upstream velocity profile for the 5L domain length, right: upstream velocity profile for the 10L domain length. Dotted blue: WMLES results, solid black line: analytical reference results.

The tests revealed the need to extend the computational domain in the streamwise direction up to a 10L length to allow the boundary layer to recover.

Nevertheless, according to the obtained results, we could conclude that the modelling strategy could reproduce the real physics provided that a sufficient streamwise domain length is prescribed. The required length will depend on the urban geometry and the degree of distortion that it introduces to the ABL.

**Evaluation**

Several visualizations have been developed to show the results of the simulations for both use cases. These visualizations were shared with the Civil Protection of Lisbon City Council.

Several postprocessing techniques have been developed to easily visualize the raw data generated by the model.

Those techniques are intended to help in the elaboration of emergency plans by the civil protection services and for sensor and other urban control systems deployment planning. To do so, animations of the temporal evolution of the hazardous substances are presented for the weather conditions specified in the report.

On the other hand, two other features allow visualizing the dispersion cloud and the pollutant concentration levels at the pedestrian level for a given instant. The data can be opened with Google Earth, allowing interactive navigation. All these methods allow predicting the first and most affected areas in the event of an accident.

On the other hand, a proof of concept of virtual sensors deployment has also been carried out to show the capabilities of the model regarding real sensor deployment planning.

Pollutant concentration section at the pedestrian level through Google Earth

This technique allows the visualization of the pollution concentration intensity in ppm at the pedestrian level. The model results are overlapped on the urban geometry and allow easy navigation through Google Earth (GE). The data is contained in a .kmz GE file that can be downloaded through the following links for both areas of study:

Alcântara:

Roma-Areeiro:

In Figure 51 and Figure 52 are presented respectively the pollutant concentration levels in ppm at pedestrian level in Alcântara and Roma-Areeiro.

###### Figure 51. Pollutant concentration levels in ppm at pedestrian level in the Alcântara area

###### Figure 52. Pollutant concentration levels in ppm at pedestrian level in the Roma-Areeiro train station area

Pollutant concentration iso-surfaces through Goolge Earth

This methodology allows to easily visualize the affected areas which are enclosed within a given concentration threshold. The concentration iso-surfaces are generated for different concentration intensities, namely 1 ppm, 10 ppm, and 100 ppm. All the volume enclosed within these surfaces features is exposed to a pollutant concentration equal to or higher than the surface value.

The data is also embedded in a kmz file to easily navigate it through Google Earth. In the following links, the iso-surfaces can be downloaded for both use cases and the three concentrations:

Alcântara:

- 1 ppm: http://dados.cm-lisboa.pt/dataset/b738cc34-4375-457d-848c-d181c82fce30/resource/ddd6573b-38b4-452a-a09d-93c1bf15075f/download/alcantaraiso-1ppm.kmz
- 10 ppm: http://dados.cm-lisboa.pt/dataset/b738cc34-4375-457d-848c-d181c82fce30/resource/7fff210f-4c35-486a-a829-47fb36e5e03d/download/alcantaraiso-10ppm.kmz
- 100 ppm: http://dados.cm-lisboa.pt/dataset/b738cc34-4375-457d-848c-d181c82fce30/resource/b7614c49-1e9e-45df-8138-52bd3ede59b4/download/alcantaraiso-100ppm.kmz

Roma-Areeiro:

- 1 ppm: http://dados.cm-lisboa.pt/dataset/b738cc34-4375-457d-848c-d181c82fce30/resource/fb4b7d33-1acb-4cae-ba4a-203524cae30c/download/roma-areeiroiso-1ppm.kmz
- 10 ppm: http://dados.cm-lisboa.pt/dataset/b738cc34-4375-457d-848c-d181c82fce30/resource/7fff210f-4c35-486a-a829-47fb36e5e03d/download/alcantaraiso-10ppm.kmz
- 100 ppm: http://dados.cm-lisboa.pt/dataset/b738cc34-4375-457d-848c-d181c82fce30/resource/b7614c49-1e9e-45df-8138-52bd3ede59b4/download/alcantaraiso-100ppm.kmz

The iso-surfaces for the three pollutant concentration levels are displayed for Alcântara in Figure 53 (1 ppm), Figure 54 (10 ppm), and Figure 55 (100 ppm).

###### Figure 53. Iso-surface of 1 ppm pollutant concentration level in Alcântara following a natural gas leakage

###### Figure 54. Iso-surface of 10 ppm pollutant concentration level in Alcântara following a natural gas leakage

###### Figure 55. Iso-surface of 100 ppm pollutant concentration level in Alcântara following a natural gas leakage

The pollutant dispersion cloud for the use case after a train accident in Roma-Areeiro was also computed.

The external cloud surface is an iso-surface concentration, considering three concentration levels: 1) iso 1 ppm (Figure 56); iso 10 ppm (Figure 57); and iso 100 ppm (Figure 58).

That means that for the iso-surface of 100 ppm, all the enclosed volume has a pollutant concentration of 100 ppm or more.

###### Figure 56. Iso-surface of 1 ppm pollutant concentration level in Roma-Areeiro train station area following a nitrobenzene spill caused by a cargo train accident

###### Figure 57. Iso-surface of 10 ppm pollutant concentration level in Roma-Areeiro train station area following a nitrobenzene spill caused by a cargo train accident

###### Figure 58. Iso-surface of 100 ppm pollutant concentration level in Roma-Areeiro train station area following a nitrobenzene spill caused by a cargo train accident

Pollutant dispersion animation through pollutant concentration

The following animations are intended to easily visualize the temporal evolution of the pollution clouds in the use case surrounding areas.

Two different animations have been generated, a top view of the urban areas of Roma Areeiro and Alcântara in which the pollutant concentration evolution is shown at street level.

In the second set of animations, the temporal evolution of the pollutant cloud coloured by concentration levels in PPM is displayed through an isometric view.

Alcântara

Simulation of the temporal evolution (in seconds) of Natural Gas concentration levels in PPM at pedestrian level following a gas leakage from the sanitation network at Alcântara area. The wind direction is assumed to be East, the most frequent wind direction in the area. The simulation methodology is based on wall-modeled Large Eddy Simulation implemented in the Barcelona Supercomputing Center’s finite-element code Alya. For this particular case, Vrema’n subgrid model together with an equilibrium Atmospheric Boundary Layer wall function has been used (Owen et al., 2020). Regarding numerical discretization, low dissipation schemes for finite elements have been used (Lehmkuhl, Houzeaux, et al., 2019). The computational mesh features 140 million elements and the simulation has been run in MareNostrum IV supercomputer using 2016 CPU, totaling 500.000 CPU-hours (see http://lisboaaberta.cm-lisboa.pt/index.php/pt/apps-e-analitica/ucd-lab for details).

- Dispersion cloud: https://www.youtube.com/watch?v=tz4G52Kiu4g

Simulation of the temporal evolution (in seconds) of Natural Gas dispersion cloud following a gas leakage from the sanitation network at Alcântara area. The cloud coloring corresponds to the concentration levels in PPM. The wind direction is assumed to be East, the most frequent wind direction in the area. The simulation methodology is based on wall-modeled Large Eddy Simulation implemented in the Barcelona Supercomputing Center’s finite-element code Alya. For this particular case, Vrema’n subgrid model together with an equilibrium Atmospheric Boundary Layer wall function has been used (Owen et al., 2020). Regarding numerical discretization, low dissipation schemes for finite elements have been used (Lehmkuhl, Houzeaux, et al., 2019). The computational mesh features 140 million elements and the simulation has been run in MareNostrum IV supercomputer using 2016 CPU, totaling 500.000 CPU-hours (see http://lisboaaberta.cm-lisboa.pt/index.php/pt/apps-e-analitica/ucd-lab for details).

Roma-Areeiro

Simulation of the temporal evolution (in seconds) of nitrobenzene concentration levels in PPM at pedestrian level following a train crash at Roma-Areeiro train station. The wind direction is assumed to be East, the most frequent wind direction in the area. The simulation methodology is based on wall-modeled Large Eddy Simulation implemented in the Barcelona Supercomputing Center’s finite-element code Alya. For this particular case, Vrema’n subgrid model together with an equilibrium Atmospheric Boundary Layer wall function has been used (Owen et al., 2020). Regarding numerical discretization, low dissipation schemes for finite elements have been used (Lehmkuhl, Houzeaux, et al., 2019). The computational mesh features 120 million elements and the simulation has been run in MareNostrum IV supercomputer using 2016 CPU, totalling 200.000 CPU-hours (see http://lisboaaberta.cm-lisboa.pt/index.php/pt/apps-e-analitica/ucd-lab for details).

- Dispersion cloud

Simulation of the temporal evolution (in seconds) of Nitrobenzene dispersion cloud following a train crash at Roma-Areeiro train station. The cloud colouring corresponds to the concentration levels in PPM. The wind direction is assumed to be East, the most frequent wind direction in the area. The simulation methodology is based on wall-modelled Large Eddy Simulation implemented in the Barcelona Supercomputing Center’s finite-element code Alya. For this particular case, Vrema’n subgrid model together with an equilibrium Atmospheric Boundary Layer wall function has been used (Owen et al., 2020). Regarding numerical discretization, low dissipation schemes for finite elements have been used (Lehmkuhl, Houzeaux, et al., 2019). The computational mesh features 120 million elements and the simulation has been run in MareNostrum IV supercomputer using 2016 CPU, totalling 200.000 CPU-hours (see http://lisboaaberta.cm-lisboa.pt/index.php/pt/apps-e-analitica/ucd-lab for details).

Simulation of the temporal evolution (in seconds) of Nitrobenzene dispersion cloud following a train crash at Roma-Areeiro train station. The cloud colouring corresponds to the concentration levels in PPM. The wind direction is assumed to be East, the most frequent wind direction in the area. The simulation methodology is based on wall-modelled Large Eddy Simulation implemented in the Barcelona Supercomputing Center’s finite-element code Alya. For this particular case, Vrema’n subgrid model together with an equilibrium Atmospheric Boundary Layer wall function has been used (Owen et al., 2020). Regarding numerical discretization, low dissipation schemes for finite elements have been used (Lehmkuhl, Houzeaux, et al., 2019). The computational mesh features 120 million elements and the simulation has been run in MareNostrum IV supercomputer using 2016 CPU, totalling 200.000 CPU-hours (see http://lisboaaberta.cm-lisboa.pt/index.php/pt/apps-e-analitica/ucd-lab for details).

Pollutant concentration sensors

This last technique is intended to help in the planning of a sensor network deployment. Three different numerical probes have been placed within the two simulation domains to emulate pollutant sensors.

The purpose of this method is to evaluate the suitability of different sensor locations to detect dangerous pollutant concentrations.

Below, the probe location map together with the simulation concentration measurements are shown for each use case.

Alcântara

###### Figure 59. Probe locations in Alcântara

###### Figure 60. Temporal evolution of natural gas concentration levels in PPM at P1 (a), P2 (b), and P3 (c) located in Alcântara.

Roma-Areeiro

###### Figure 61. Probe locations in Roma-Areeiro

###### Figure 62. Temporal evolution of nitrobenzene concentration levels in PPM at P1 (a), P2 (b), and P3 (c) located in Roma-Areeiro

**References**

Calafell, J., Trias, F. X., Lehmkuhl, O., & Oliva, A. (2019). A time-average filtering technique to improve the efficiency of two-layer wall models for large eddy simulation in complex geometries. Computers & Fluids, 188, 44–59. https://doi.org/https://doi.org/10.1016/j.compfluid.2019.03.026

García-Sánchez, C., van Beeck, J., & Gorlé, C. (2018). Predictive large eddy simulations for urban flows: Challenges and opportunities. Building and Environment, 139, 146–156. https://doi.org/https://doi.org/10.1016/j.buildenv.2018.05.007

Hine, T. B. (1924). The Rate of Evaporation of Liquids in a Current of Air. Phys. Rev., 24(1), 79–91. https://doi.org/10.1103/PhysRev.24.79

Kurppa, M., Hellsten, A., Auvinen, M., Raasch, S., Vesala, T., & Järvi, L. (2018). Ventilation and Air Quality in City Blocks Using Large-Eddy Simulation—Urban Planning Perspective. Atmosphere, 9(2). https://doi.org/10.3390/atmos9020065

Lehmkuhl, O., Houzeaux, G., Owen, H., Chrysokentis, G., & Rodriguez, I. (2019). A low-dissipation finite element scheme for scale resolving simulations of turbulent flows. Journal of Computational Physics, 390, 51–65. https://doi.org/https://doi.org/10.1016/j.jcp.2019.04.004

Lehmkuhl, O., Piomelli, U., & Houzeaux, G. (2019). On the extension of the integral length-scale approximation model to complex geometries. International Journal of Heat and Fluid Flow, 78, 108422. https://doi.org/https://doi.org/10.1016/j.ijheatfluidflow.2019.108422

Li, X.-X., Britter, R. E., Koh, T. Y., Norford, L. K., Liu, C.-H., Entekhabi, D., & Leung, D. Y. C. (2010). Large-Eddy Simulation of Flow and Pollutant Transport in Urban Street Canyons with Ground Heating. Boundary-Layer Meteorology, 137(2), 187–204. https://doi.org/10.1007/s10546-010-9534-8

Liu, Y. S., Miao, S. G., Zhang, C. L., Cui, G. X., & Zhang, Z. S. (2012). Study on micro-atmospheric environment by coupling large eddy simulation with mesoscale model. Journal of Wind Engineering and Industrial Aerodynamics, 107–108, 106–117. https://doi.org/https://doi.org/10.1016/j.jweia.2012.03.033

Lozano-Durán, A., & Jiménez, J. (2014). Effect of the computational domain on direct simulations of turbulent channels up to Re τ = 4200. Physics of Fluids, 26(1), 011702. https://doi.org/10.1063/1.4862918

Lu, L., Zhang, X., Yan, Y., Li, J., & Zhao, X. (2014). Theoretical Analysis of Natural-Gas Leakage in Urban Medium-pressure Pipelines.

Owen, H., Chrysokentis, G., Avila, M., Mira, D., Houzeaux, G., Borrell, R., … Lehmkuhl, O. (2020). Wall-modeled large-eddy simulation in a finite element framework. International Journal for Numerical Methods in Fluids, 92(1), 20–37. https://doi.org/https://doi.org/10.1002/fld.4770

Pastrana, D., Cajas, J. C., Lehmkuhl, O., Rodríguez, I., & Houzeaux, G. (2018). Large-eddy simulations of the vortex-induced vibration of a low mass ratio two-degree-of-freedom circular cylinder at subcritical Reynolds numbers. Computers & Fluids, 173, 118–132. https://doi.org/https://doi.org/10.1016/j.compfluid.2018.03.016

Rodriguez, I., Lehmkuhl, O., Soria, M., Gómez, S., Domínguez-Pumar, M., & Kowalski, L. (2019). Fluid dynamics and heat transfer in the wake of a sphere. International Journal of Heat and Fluid Flow, 76, 141–153. https://doi.org/https://doi.org/10.1016/j.ijheatfluidflow.2019.02.004

Tominaga, Y., Mochida, A., Yoshie, R., Kataoka, H., Nozu, T., Yoshikawa, M., & Shirasawa, T. (2008). AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings. Journal of Wind Engineering and Industrial Aerodynamics, 96(10), 1749–1761. https://doi.org/https://doi.org/10.1016/j.jweia.2008.02.058

Toparlar, Y., Blocken, B., Maiheu, B., & van Heijst, G. J. F. (2017). A review on the CFD analysis of urban microclimate. Renewable and Sustainable Energy Reviews, 80, 1613–1640. https://doi.org/https://doi.org/10.1016/j.rser.2017.05.248

Vázquez, M., Houzeaux, G., Koric, S., Artigues, A., Aguado-Sierra, J., Arís, R., … Valero, M. (2016). Alya: Multiphysics engineering simulation toward exascale. Journal of Computational Science, 14, 15–27. https://doi.org/https://doi.org/10.1016/j.jocs.2015.12.007

Yaghoobian, N., Kleissl, J., Tha Paw, K. U., Yaghoobian, N., Kleissl, J., & Paw U, K. T. (2014). An Improved Three-Dimensional Simulation of the Diurnally Varying Street-Canyon Flow. 153, 251–276. https://doi.org/10.1007/s10546-014-9940-4