2.1 Model configuration

All runs are forced by ERA5 reanalysis (^{Hersbach et al., 2017}) at spatial and temporal resolutions of 30 km and 3 h, respectively. The control simulation (c) uses the following schemes: Mellor-Yamada-Janjic for PBL, Eta similarity for surface layer, Noah for land surface, Thompson for microphysics, RRTMG for radiation, and no urban canopy. Other configurations were selected by changing one scheme at a time as shown in Table I and Figure S1 of the supplementary material. Included PBL schemes are: p¹: Yonsei University (YSU, ^{Hong et al., 2006}), p²: Mellor- Yamada-Janjic (MYJ, ^{Janjić, 1994}), p⁶: Mellor-Yamada ^{Nakanishi Niino (MYYN, Nakanishi and Niino, 2006}), p⁷: Asymmetric Convection Model 2 Scheme (ACM2, ^{Pleim, 2007}), and p⁸: Bougeault-Lacarrere Scheme (^{Bougeault and Lacarrere, 1989}). Included surface layer schemes: s¹: Revised MM5 (^{Jiménez et al., 2012}), s²: Eta Similarity (^{Janjić, 1994}), s⁵: MYNN and s⁷: Pleim-Xiu (^{Pleim, 2006}). Included land surface schemes are: l²: Unified Noah Land Surface Model (Noah, ^{Tewari et al., 2004}) and l⁴: Noah-MP Land Surface Model (Noah-Mp, ^{Niu et al., 2011}). Included radiation schemes are: r⁴: RRTMG (^{Iacono et al., 2008}), r⁷: Fu-Liou-Gu (^{Gu et al., 2011}), and r¹⁴: RRTMG-K (^{Baek, 2017}). Included microphysics schemes: m¹⁰: Morrison 2-moment (^{Morrison et al., 2009}), and m⁵²: P3 (^{Morrison and Milbrandt, 2015}). Included urban schemes: u¹: Urban Canopy Model (^{Chen et al., 2011}), u²: BEP (^{Martilli et al., 2002}), and u³: BEP + BEM (^{Salamanca et al., 2010}).

As the MABA is a highly urbanized area, an improvement in model performance could be expected with the activation of an urban scheme. For that reason, PBL schemes that are compatible with urban options (in WRF 4.2.1: Mellor-Yamada-Janjic and BouLac) are tested. Although Yonsei University (YSU) and Asymmetric Convection Model 2 (ACM2) are not compatible with the Building Effect Parameterization (BEP) and the Building Energy Model (BEP + BEM) urban schemes in the WRF version used in this work, these are still included as they are widely used in similar studies for air quality applications (e.g., ^{Cuchiara et al., 2014}; ^{Banks and Baldasano, 2016}; ^{Mohan and Gupta, 2018}). BEP (^{Salamanca and Martilli, 2010}) is a 3D urban canopy model which accounts for multiple building parameters allowing for higher constructions than the first model level. BEP + BEM (^{Martilli et al., 2002}) is a multilayer building energy model composed of a building energy model coupled to BEP. BEM accounts for the impacts of anthropogenic heat emissions on the urban environment and estimates the cooling/heating energy demand due to air-conditioning systems. When a UCM is not used, the Noah and Noah-MP LSMs use the Bulk scheme to represent the urban surface in the WRF model. Although the PBL, urban, and land surface model schemes are expected to have a greater impact on model results, we analyze the sensitivity to other physics options (such as radiation and microphysics) as this is the first exploratory study of this kind in the area and these can also affect model performance as found in other works (e.g., ^{Borge et al., 2008}).

2.2 Domain

The MABA is composed of the city of Buenos Aires and 24 districts of greater Buenos Aires. It has a population of around 15 million inhabitants over 3800 km² (^{INDEC, 2010}) of flat terrain and it is surrounded by non-urban areas and the La Plata River on its east side. It is the third most populated megacity in Latin America (^{UN, 2019}). Three modeling nested domains with resolutions of 15 km (111 × 101 grid points), 3 km (121 × 101 grid points), and 1 km (133 × 127 grid points) (Fig. 1) are used with the outermost domain (d01) comprising the whole Buenos Aires province and the innermost domain (d03) covering the whole MABA. Eighty terrain-following hybrid vertical levels near the ground are used. Twelve of them cover the lowest 3 km from the surface with the first level at 24 m. Land uses for this region come from satellite data (Moderate Resolution Imaging Spectroradiometer [MODIS] at a 30” horizontal resolution). Unfortunately, there is only one urban class available for the area at the moment of this study.

Fig. 1 WRF modeling nested domains (left) and map of the MABA (right) with the meteorological (red crosses and blue crosses) and air quality (yellow crosses) stations. Source: Google Maps. 

2.3 Period

One week in winter and one week in spring of 2012 were chosen to cover different meteorological conditions with enough available data on air pollutants concentration for future air quality model performance validation studies. The first two days of each period are considered as a spin-up. The winter week (July 28 to August 4) started with a high-pressure anomaly on the Atlantic Ocean and a low-pressure anomaly on the Pacific Ocean. This led to a rotation of wind direction and the week ended with precipitation. On the other hand, on the spring period (November 10 to 17) there was a high-pressure disturbance that lasted the whole week with mostly southeastern winds and dry air over the city, as discussed in ^{Luque et al. (2021)}. Synoptic maps are provided in the supplementary material. All dates presented in this work are expressed in UTC time while the MABA has a local time of UTC-3.

2.4. Model performance evaluation

To assess model performance, four meteorological variables that are relevant for air quality purposes are considered (^{Ziomas et al., 1995}; ^{Harkey et al., 2015}): 2-m temperature (T) and water vapor mixing ratio (Qv), and 10-m wind speed (Ws) and direction (Wd). Modeled values of T, Qv, Ws, and Wd are statistically compared with those observed in the Aeroparque meteorological station (AEP, see Fig. 1). This station is the closest one to the three available air quality stations whose data will be used for validation of future air quality simulations and the most representative one. The following statistics are computed: bias, mean absolute error (Mae), root mean square error (Rmse), Correlation (Corr), and index of agreement (Indexagr2) (^{Emery et al., 2001}). This last metric is calculated as the ratio between the Rmse and the sum of the difference between each prediction and the observed mean, and the difference between each observation and the observed mean. It is a measure of the match between the departure of each prediction and each observation from the observed mean. Configurations with the lowest errors and higher Corrs and Indexagr2 are selected to be further analyzed.

For these configurations, the relationship between errors in Wd and Ws is studied following the methodology by ^{Jiménez and Dudhia (2013)}. In addition, the planetary boundary layer height (PBLH) is explored. PBLH results are contrasted with observations from sounding data available only at 12:00 UTC taken in the Ezeiza Meteorological station (EZE, see Fig. 1). In order to compare PBLH results obtained from the different configurations, modeled and observed values are recalculated following the methodology described in ^{Nielsen-Gammon et al. (2008)} and used in other works (e.g., ^{Miao et al., 2022}; ^{Yan et al., 2022}). This method estimates the PBLH as the height where a “critical inversion occurs”. This height is defined in ^{Marsik et al. (1995)} as the level where potential temperature is 1.5 K higher than the minimum value it presents in the PBL.

In order to assess model performance throughout the whole MABA domain, performance metrics are then computed at other meteorological stations for the best performing configurations, which are presented in the supplementary material.

3.1 Hourly evolution

Figure 2 presents the hourly evolution of modeled and observed values of Ws, Wd, T, and Qv. Observed Ws values vary between 0 and 10 m s^-1, with the lowest ones taking place during the winter week (Fig. 2a, b). Temporal variations for this variable are different between each week and most configurations follow these patterns except for those with complex urban schemes (u², u³, and p⁸u³) that show lower wind intensities. During winter week winds come mostly from the SE and SW except for July 31 and August 2, when they come from the NE and N, with very low intensities (Fig. 2c). During the spring week, winds come mostly from the SE (Fig. 2d). This variable is also quite well reproduced by most configurations during both weeks.

Fig. 2 Temporal evolution of modeled and observed wind speed (Ws [m s^-1]), direction (Wd [º]), temperature (T [K]) and vapor mixing ratio [Qv, g kg^-1] during winter (left column) and spring (right column) weeks. Hours are written in UTC and the MABA has a local hour of UTC-3. (a) and (b) Ws, (c) and (d) Wd, (e) and (f) T, (g), and (h) Qv. 

Observed T values (Fig. 2e, f) display different behaviors between the two periods, with a clearer diurnal cycle during spring. Most configurations follow observations, even though for some of them differences tend to become larger over time showing the expected decrease in model skill as time of forecast increases (^{Buizza and Leutbecher, 2015}).

During the winter week, observed values of Qv rapidly grow until August 1 (Fig. 2g), when precipitation starts reaching its maximum value. This behavior is due to the fact that on July 31 there is a high-pressure zone in front of the coastline of Uruguay centered around 35º S and 55º W (see Figs. S2 to S7 in the supplementary material). The circulation associated with the high-pressure anomaly advects wetter air towards the Buenos Aires area. After that, variations are smaller. During the spring week, an increase in Qv is also seen (Fig. 2h). During both periods, all simulations reproduce well the time evolution of Qv but underestimate its values between -16 and -5% depending on the configuration, mainly in spring.

3.2 Performance metrics

In order to select the best performing configurations, Figure 3 presents the bias, Mae, Rmse, Corr and Indexagr2 computed for each variable, week and configuration. For values of each statistic see section S3 in the supplementary material. For Ws, configurations present similar performance metrics during both weeks (Fig. 3a), but relative error values are smaller in spring (not shown). Relatively large values of Corr (> 0.77) and Indexagr2 (> 0.88) are obtained for most schemes. Configurations with complex urban schemes (which include the label u² or u³: u², u³, p⁸u², and p⁸u³) present the worst performance metrics with Mae 20-30% higher than the other configurations and lower Correlations, specially u² and u³ during the spring week, which present Correlation values below 0.3. While a better performance with these configurations could be expected in the highly urbanized MABA, some works show that complex urban schemes behave better when local urban parameters suited for the particular city under study are used (e.g., ^{Kuik et al., 2016}; ^{Shen et al., 2019}). In this first stage of the model implementation, default WRF urban parameters are used. A sensitivity test to the multiple parameters (e.g., building height and albedo, among others) used to configure the urban schemes is beyond the scope of this work due to computational constraints and the lack of a set of available reference values for these parameters in MABA. The larger underestimation for this variable obtained with configurations with BEP and BEP + BEM (u², u³, p⁸u², and p⁸u³) could be due to an overestimation of the urban fraction in AEP as a consequence of using default urban parameters. The other configurations without complex urban schemes present bias, Mae and Rmse in the ranges of -1.6 to 0.5 m s^-1 (with relative values of -28.2 to 1.3%), 0.9 to 1.8 m s^-1 (17.6 to 33.2%), and 1.5 to 2.1 m s^-1 (27.7 to 53.2%), respectively. During both weeks, configurations with a simple urban scheme (those labelled with u¹, p⁸u¹ and l⁴u¹) have the lowest errors. In the spring, p⁸u¹ has the best performance as it also presents the higher Corr and Indexagr2.

Fig. 3 Statistics (vertical axis) for Ws (m s^-1), Wd (º), T (K) and Qv (g kg^-1) for both weeks and all configurations (horizontal axis). Index of agreement and Correlation in green and errors in red, white and blue. (a) Ws, (b) Wd, (c) T, and (d) Qv. 

Wd (Fig. 3b) is counter-clockwise biased by almost all configurations during both periods with bias ranging between -23.4º and -4.2º (-16.6 and -2.9% with stronger underestimations during spring), Mae between 18.7º-52.2º (17.0-37.0%) and Rmse in the range 40.9º-75º (37.2-53.4%). During the winter week, Mae and Rmse values are relatively higher. The configuration with the MYNN PBL scheme (p⁶) has a bias close to 0% but this is clearly due to a compensation of errors in the calculation, as this configuration also has the highest Rmse and Mae. Both Corr and Indexagr2 are higher during the spring with values greater than 0.6. Configurations with the Noah-MP land surface scheme (l⁴, l⁴m¹⁰ and l⁴m⁵²), with the exception of the combination with the simple urban scheme (l⁴u¹), have clearly the lowest errors during both weeks with similar Corr and Indexagr2 to those of other configurations.

For T (Fig. 3c), Corr and Indexagr2 are higher during the spring week than during the winter week with values above 0.9. This is expected as a better representation of the diurnal cycle is observed in the spring period (Fig. 2e, f). Still, all Corr and Indexagr2 have values over 0.7 in the winter week. For both periods, bias, Mae and Rmse values range from -1.4 to 1.2 K (-0.5-0.4%), 0.5 K to 1.5 K (0.2-0.5%) and 1.8K to 3.4 K (0.6-1.2%) respectively. Configurations with Noah-Mp land surface scheme without urban scheme (l⁴, l⁴m¹⁰ and l⁴m⁵²) slightly underestimate T while the only configuration with Noah-Mp and urban scheme (l⁴u¹) overestimates it. This is most likely due to the activation of the urban scheme. Other configurations share this characteristic and do not exhibit a similar positive bias, but they also have a different land surface model (u¹) and a different PBL scheme (p⁸u¹). The use of urban schemes improves the performance during the spring week, but not during the winter one. A possible explanation for this is that the impact of urban environments is larger during the nocturnal period, when the impact of accumulated heat from buildings and paved surfaces is more noticeable (^{Argüeso et al., 2014}). Therefore, it is expected that the impact of the presence of urban structures would be more noticeable during the warmer season, when the solar radiation is also more abundant.

Qv is underestimated during both weeks, as observed in Figure 2g, h with bias values in the range of -1.1 and -0.3 g kg^-1 (-4.65 and -1.39%) (Fig. 3d). Mae and Rmse are both in the range of 1.3-2.5 g kg^-1. Corr and Indexagr2 are higher in the winter week than in the spring one with values over 0.9. During the spring period, on the other hand, Corr and Indexagr2 values are close to 0.7 with the exception of configurations with Noah-Mp (l⁴, l⁴m¹⁰ and l⁴m⁵²), which have similar values to those obtained in the winter week, and p⁷, that presents lower values than other configurations.

In order to gain preliminary insight on the impact of each physical scheme on surface meteorological variables, configurations that share all characteristics except one are compared: c (control run), m¹⁰, and m⁵² differ only in their microphysics scheme (Thompson, Morrison 2-moment, and P3, respectively) and as can be seen in Figure 3, the four configurations present similar behavior, suggesting that this scheme does not significantly affect the model performance. The same is observed with the radiation scheme, as configurations c, r⁷ and r¹⁴ differ only in their microphysics scheme (RRTMG, Fu-Liu-Gu, and RRTMG-K, respectively) and they present similar results.

Configurations c, u¹, u² and u³ differ only in their urban scheme (Bulk scheme, simple layer, BEP, and BEP + BEM, respectively). Contrary to what happens with the microphysics and radiation parametrizations, here we see a clear effect. Results of the three configurations with urban schemes (u¹, u², and u³) change mostly in the case of Ws, giving lower values than the control configuration due to the increase in surface roughness length and blocking effects associated to buildings. They also affect results for T. Contrary as expected, the activation of the urban scheme lowers its values during both weeks.

Configurations c, p¹, p⁶, p⁷, and p⁸ differ on their PBL and surface scheme (MYJ + Eta similarity, MYNN, ACM2+Pleim, Xiu and BouLac + Revised MM5, respectively). Control simulation differs from the others with higher Ws and T values, but these subsets of simulations behave similarly. The exception is p⁷, which has a better performance for Wd during the winter week, but a lower one for Qv and T during the spring week.

Lastly, configurations c and l⁴ differ only in their land surface model (Noah for c and Noah MP for l⁴). Therefore, minimum differences are expected between the results obtained with both configurations at an urban grid cell such as AEP; however, l⁴ shows a slight improvement in the representation of Qv and Wd compared with c, especially during the spring week. A possible explanation for this is that since both variables depend on their values in other cells, an impact from remote non-urban cells could be affecting the results found in AEP.

These comparisons suggest that the urban scheme has a greater influence on the model performance while radiation and microphysics schemes have a negligible effect. However, since the response to different physical schemes is not linear, more configurations should be compared with each other at more weather stations to have a clearer conclusion regarding this matter.

Almost all configurations have good performances for T and Qv. This is consistent with previous works where the dynamics of these variables are easier to represent by the model than Ws and Wd (e.g., ^{Cuchiara et al., 2014}; ^{Banks and Baldasano, 2016}; ^{Mohan and Gupta, 2018}). For this reason, the selection of configurations for further inspection is based on model performance to estimate the wind variables.

It should be noted that while Wd is mostly influenced by synoptic/mesoscale conditions, Ws can also be influenced by local dynamics that depend on model configuration. This could explain that while dynamics of the predominant synoptic/mesoscale features might not be too different between simulations, the wind intensity might be. The best performance to estimate Ws and l⁴ is given by l⁴m¹⁰, while l⁴m⁵² presents the best performance for Wd. Configurations p¹ and l⁴u¹ are also selected to analyze the impact of different PBL schemes (a local scheme and a non-local one, respectively) and whether the activation of a simple urban scheme plays a role in performance.

Figure 4 shows the differences between hourly observed and modeled values for all variables in each of the six configurations. For Ws, most differences are within the range -2-2 m s^-1 and behave similarly among configurations (Fig. 4a, b). During both weeks, p⁸u¹ usually presents slightly higher Ws values than the rest (typically by 1 m s^-1, which represents around 25 or 17% of the observed mean, depending on the week). This leads to a somewhat larger overestimation of minimums but a smaller underestimation of maximums. On the other hand, p¹ (configured with YSU as PBL scheme) presents lower Ws values (around 1 m s^-1); therefore, it overestimates less minimums values, but underestimates more maximums values. Configurations with Noah-Mp are an intermediate case. During some periods of a few hours, the difference between model results and observations are larger than 2 m s^-1, which represents around 34 or 50% of the observed mean depending on the week. For all configurations, the largest underestimation (6 m s^-1) coincides with the global observed maximum of Ws (see Fig. 2a).

Fig. 4 Differences between modeled values of Ws (m s^-1), Wd (º), T (K) and Qv (g kg^-1) with observations for chosen configurations for the winter (left column) and spring week (right column). The range in which most differences lie is colored in grey. (a) and (b) Ws, (c) and (d) Wd, (e) and (f) T, (g) and (h) Qv. 

Differences between modeled and observed Wd values for the different configurations are also similar to each other and within the range -50-50º (Fig. 4c, d). Larger differences are observed around August 2 at 00:00 UTC when Ws values are close to 0 m s^-1 and the midday of November 13 and 16, when configurations with Noah-MP predict Ws values lower than 2 m s^-1.

In the case of T (Fig. 4e, f), most differences with observations are in the range -2 to 2 K but larger differences among the selected configurations are observed, compared to those of wind variables. During both weeks, l⁴u¹ shows larger (and mainly positive) differences between modeled and observed T values, while best performing schemes for Wd (l⁴, l⁴m¹⁰, and l⁴m⁵²) present lower (and mostly negative) differences, as also shown in Figure 4c. During the spring week, l⁴u¹ shows a larger overestimation of temperature, which increases with time during the last three days of the period. On the other hand, configurations p¹ and p⁸u¹ usually present smaller differences. To understand why configuration l⁴u¹ (Noah-Mp + single layer urban scheme) overestimates observations and presents values between 1 and 4 K higher than the others, surface heat fluxes are explored (see section S4 in the supplementary material). A higher sensible heat flux for this configuration (l⁴u¹) is expected, but it is found that it presents values similar to the other, all with higher values during the spring week. To understand why it still presents higher T, surface temperature is analyzed, since sensible heat flux depends on the temperature difference between the surface and the air above it. If surface temperature is higher in this configuration, it could explain why it presents similar sensible heat values as the other configurations, but with a higher T. This is the case, and is most likely due to the activation of the urban scheme.

For Qv (Fig. 4g, h), differences between model results and observations are mostly negative and vary similarly. During the winter week, most differences range between -1.5 and 1.0 g kg^-1. Configurations p¹ and p⁸u¹ present larger differences with observations, especially during the spring week where a greater underestimation (up to 6 g kg^-1) is observed. From Figure 2h, it can be seen that in this week the largest underestimations occur around minimums.

To better understand the relationship between wind direction errors and wind speed, the root mean square error of Wd is calculated for different subsets of data, eliminating those with observed Ws values below a given threshold, which varies between 0 and 5 m s^-1 covering most of the dataset in order to find if a change in behavior exists at some point and how configurations behave after it. This is done for the six selected configurations (Fig. 5). Rmse for all configurations in both weeks is larger for lower Ws threshold values, consistent with results by ^{Jiménez and Dudhia (2013}). During the winter week (Fig. 5a), the Wd Rmse varies from 44º to 18º for Noah-Mp configurations, 73º to 20º for p¹, and 88º to 24º for p⁸u¹. A very small variation of Wd Rmse is observed for Ws threshold values above 2.6 m s^-1. The most abrupt decrease (70%) is presented by p⁸u¹, highlighting its large sensitivity to this variable for low Ws values. Note that around 35% of wind speed measurements during this week have wind speed values below 2.6 m s^-1. In turn, during the spring week (Fig. 5b), the Wd Rmse varies form 38º to 19º, representing a change of around 40%) for all configurations. This change in Rmse occurs when the threshold is around 3.1 m s^-1. In this case, only 19% of the observed hourly Ws values are below this value. The larger presence of lighter winds during the winter week may explain the relatively worse results to simulate wind variables during this period.

Fig. 5 Wd (º) Rmse (left side scale) for different Ws (m s^-1) threshold and Ws frequency (m s^-1) (right side scale) for different wind speed values. Filled light grey bars show the normalized probability density function for observed Ws. (a) Winter week, (b) spring week. 

3.3. Planetary boundary layer height

Figure 6 shows modeled PBLH values at AEP and EZE stations (Fig. 1) during both weeks. Only five sounding data points are available for EZE (black dots) during the analyzed period, and they are included only for comparison in the EZE plots. For each week, the temporal variations of the PBLH are similar at the two stations, with maximum diurnal values being somewhat higher at AEP. At these sites, modeled PBLH values vary from 100 to 1500 m during the winter week and from 250 to 2000 m in the spring week. These values are a bit higher than the ones studied in ^{Mazzeo and Gassmann (1990)}, who found average maximum values of 1100 m during August and 1500 m during November. In general, configurations present similar temporal variations with higher PBLH values during the day and lower ones at nighttime, as expected. However, large differences are observed on some days. For example, on July 31, p⁸u¹ presents larger values than other configurations at the two stations. This is consistent with the fact that this configuration also presents higher wind speed values leading to more mechanical turbulence during that day. On November 16 and 17, p⁸u¹ and p¹ simulate PBLH peaks that are considerably higher (by 1000 m) than those estimated by other configurations. Although on November 17 p¹ and p⁸u¹ present stronger wind intensities that could partially explain this, no significant differences in wind speed between configurations are observed on November 16.

Fig. 6 Planetary boundary layer height for chosen configurations and observations (black dots) for winter (left column) and spring weeks (right column). (a) and (b) Aeroparque, (c) and (d) Ezeiza. 

3.4. Performance evaluation at other meteorological stations

Model performances metrics (bias, Mae, Rmse, Corr and Indexagr2) are computed for best performing configurations for Ws and Wd at AEP (p⁸u1 and l⁴, respectively) at other five meteorological stations of the MABA (San Miguel [SMG], Ezeiza [EZE], Moron [MOR], Observatorio de Buenos Aires [OBS] and Palomar [PAL]; see Fig. 1). Results are included in section S5 in the supplementary material. Both configurations present a good performance to estimate T with Corr > 0.7 and errors under 15% for all stations. Similar results are obtained for Qv with the exception of San Miguel and Moron stations during the spring week. On the other hand, model performance to estimate Ws and Wd varies with the station and the week. Ws presents an acceptable performance (Corr > 0.6 and Mae < 35%) in Ezeiza and Palomar stations during both weeks and in Moron during the spring period, but a poor one in San Miguel and Palomar during the winter week. Wd has an acceptable performance in Ezeiza and in Moron during the winter period. The worst model performance for both wind variables is obtained at OBS, which may be expected as this station is surrounded by buildings and trees that affect wind measurements (^{Mazzeo and Gassmann, 1990}).