Assessing Spatial Variation of PBL Height and Aerosol Layer Aloft in S ã o Paulo Megacity Using Simultaneously Two Lidar during Winter 2019

: This work presents the use of two elastic lidar systems to assess the horizontal variation of the PBL height (PBLH) and aerosol layer aloft in the S ã o Paulo Megacity. These two lidars performed simultaneous measurements 10.7 km apart in a highly urbanized and relatively ﬂat area of S ã o Paulo for two winter months of 2019. The results showed that the PBLH differences display diurnal variation that depends on the PBL during daytime growth phases. Cloud and sea breeze effects control most of PBLH variation. In the absence of cloud and sea breeze, the maximum difference (~300 m) occurs in the rapid development stage and is due to topographic effects. When the PBL approaches its maximum daily value, it tends to level off with respect to the topography. In addition, it was presented a method that combines elastic lidar (to detect an aerosol layer) and satellite data (to classify such a layer from Aerosol Optical Depth (AOD) and Aerosol Index (AI) information) for the detection of biomass burning events. This methodology demonstrated that the variations caused by Biomass Burning in AOD and AI enable both the detection of aerosol plumes originating from biomass burning and the identiﬁcation of their origin.

In the case of elastic lidars [13], the PBLH detection methods are based on the abrupt reduction in the backscattered signal intensity that occurs in the transition layer between the PBL and Free Atmosphere (FA). Unfortunately, such a drop in the backscattered signal intensity occurs in ideal conditions, observed during the convective period in the absence of low clouds and decoupled aerosol layers. As ideal conditions are rarely observed, several

Metropolitan São Paulo Lidar 1 (MSP1) System
The Metropolitan São Paulo Lidar 1 (MSP1) system is a coaxial ground-based multiwavelength Raman lidar located at the Nuclear and Energy Research Institute (IPEN) in the suburban region of São Paulo city (23°34′ S, 46°43′ W, 782 m asl) (Figure 1b). MSP1 operates with a pulsed Nd:YAG laser pointed towards the zenith direction, emitting radiation at 355 nm, 532 nm, and 1064 nm with a repetition frequency of 10 Hz. This system detects three Raman-shifted channels (387, 408, and 530 nm) and three elastic channels (355, 532, and 1064 nm), reaching the full overlap at 300 m above ground level [7]. During the field campaign, MSP1 was run with a temporal and spatial resolution of 1 min and 7.5 m, respectively, from 09 to 18 LT.

Metropolitan São Paulo Lidar (MSP) 2 System
The Metropolitan São Paulo Lidar 2 (MSP2) is a mobile biaxial ground-based multiwavelength Raman lidar system. MSP2 operates with a pulsed Nd: YAG laser, which emits radiation at 355 nm and 532 nm in the zenith direction and detects one elastic channel (532 nm) and one Raman-shifted channel (387 nm). Such a system reaches the full overlap at 180 m above ground level. During the field campaign, this system was allocated in the Department of the Treasury of the São Paulo State (SEFAZ) in the urban region of São Paulo city (23°32′ S, 46°37′ W, 739 m asl) (Figure 1b), operating continuously with a temporal and spatial resolution of 1 min and 7.5 m, respectively.

Suomi National Polar-Orbiting Partnership (Suomi NPP) Data
The Suomi National Polar-orbiting Partnership (Suomi NPP) is a weather satellite operated by the National Oceanic and Atmospheric Administration (NOAA). Launched in 2011, Suomi NPP is equipped with the Visible Infrared Imaging Radiometer Suite (VIIRS) and the Ozone Mapping and Profiling Suite Nadir-Mapper (OMPS-NM). The VIIRS is a whiskbroom scanner radiometer, which passively observes reflectance at visible and infrared wavelengths [35].
In this work, Aerosol Optical Depth (AOD) and Aerosol Index (AI) were obtained respectively from VIIRS data and normalized radiances using two-wavelength pairs at 340 and 378.5 nm of OMPS-NM. AI is a qualitative index that indicates the presence of aerosol layers, such as biomass-burning, desert-dust, and volcanic-ash plumes, monitoring the absorption in the above-mentioned wavelength pairs [36]. Both datasets are available at: https://earthdata.nasa.gov (accessed on 01 March 2021).   (Figure 1b). MSP1 operates with a pulsed Nd:YAG laser pointed towards the zenith direction, emitting radiation at 355 nm, 532 nm, and 1064 nm with a repetition frequency of 10 Hz. This system detects three Raman-shifted channels (387, 408, and 530 nm) and three elastic channels (355, 532, and 1064 nm), reaching the full overlap at 300 m above ground level [7]. During the field campaign, MSP1 was run with a temporal and spatial resolution of 1 min and 7.5 m, respectively, from 09 to 18 LT.

Metropolitan São Paulo Lidar (MSP) 2 System
The Metropolitan São Paulo Lidar 2 (MSP2) is a mobile biaxial ground-based multiwavelength Raman lidar system. MSP2 operates with a pulsed Nd: YAG laser, which emits radiation at 355 nm and 532 nm in the zenith direction and detects one elastic channel (532 nm) and one Raman-shifted channel (387 nm). Such a system reaches the full overlap at 180 m above ground level. During the field campaign, this system was allocated in the Department of the Treasury of the São Paulo State (SEFAZ) in the urban region of São Paulo city (23 • 32 S, 46 • 37 W, 739 m asl) (Figure 1b), operating continuously with a temporal and spatial resolution of 1 min and 7.5 m, respectively.

Suomi National Polar-Orbiting Partnership (Suomi NPP) Data
The Suomi National Polar-orbiting Partnership (Suomi NPP) is a weather satellite operated by the National Oceanic and Atmospheric Administration (NOAA). Launched in 2011, Suomi NPP is equipped with the Visible Infrared Imaging Radiometer Suite (VIIRS) and the Ozone Mapping and Profiling Suite Nadir-Mapper (OMPS-NM). The VIIRS is a whiskbroom scanner radiometer, which passively observes reflectance at visible and infrared wavelengths [35].
In this work, Aerosol Optical Depth (AOD) and Aerosol Index (AI) were obtained respectively from VIIRS data and normalized radiances using two-wavelength pairs at 340 and 378.5 nm of OMPS-NM. AI is a qualitative index that indicates the presence of aerosol layers, such as biomass-burning, desert-dust, and volcanic-ash plumes, monitoring the absorption in the above-mentioned wavelength pairs [36]. Both datasets are available at: https://earthdata.nasa.gov (accessed on 1 March 2021).

AERONET Sunphotometer
The AErosol RObotic NETwork (AERONET) [37] is the NASA sunphotometer global network that supplies automatic sun and sky scanning measurements. Using direct sun measurements, AERONET supplies both Aerosol Optical Depth (AOD) and the Ångström Exponent (Å), which gives the wavelength dependence of the AOD. Using multiangular and multispectral measurements of atmospheric radiance and applying a flexible inversion algorithm [38], the AERONET data can also supply several additional aerosol optical and microphysical parameters, such as size distribution, single-scattering albedo, and refractive index. The operating principle of this system is to acquire aureole and sky radiance observations using a large number of solar scattering angles through a constant aerosol profile and thus retrieve the aerosol size distribution, the phase function, and the AOD [39]. In this work, the AERONET sunphotometer data were measured in the SP-EACH station located in São Paulo city (23 • 48 S, 46 • 49 W, 754 m asl). These data were used to derive the aerosol size distribution and Angstrom Matrix during the BB event from 17 to 19 August 2019.

PBLH Detection
The PBLH was estimated from the 532 nm backscattered signal, so the algorithm is divided into two parts: raw lidar data pre-processing [40] and Wavelet Covariance Transform (WCT) algorithm [14]. The data pre-processing begins with the subtraction of the dark current signal (DC(z)) of the raw lidar signal (P) to reduce the influence of electrical noise. Then, the background radiation signal (BG) is removed to attenuate the influence of external sources. Finally, due to the attenuation of the lidar signal with the height, the result of the previous steps is multiplied by the square of the corresponding height (z), resulting in the Range Corrected Signal (RCS 532 ), as indicated in Equation (1): After the raw signal pre-processing, the second part of the algorithm is performed, where the WCT algorithm is applied. Firstly, the covariance (W(a, b)) between the average RCS 532 (z) obtained for one hour (RCS 532 (z)) and Haar function (h z−b a ), the motherwavelet is done: where z is the height above the ground, z i and z f are the lower and upper limit of the RCS 532 (z), respectively, and the respective values of dilatation and transition-related to mother-wavelet are given by a and b. Considering a previous study by Moreira et al. [41], which was held in São Paulo city, the parameters a and b received the values 200 and 40 m, respectively. Then, the height z where the maximum W(a, b) occurs is classified as PBLH:

PBLH Levelness
The Levelness number (L) was defined by Stull [42] as: where ∆z i and ∆z T are PBLH and topography differences, respectively. The ratio L provides information about how the PBLH varies horizontally with respect to the topographic variations. If L < 0, the PBLH varies in an opposite way to the topography, that is, the higher terrains have lower PBLH values and vice-versa. L = 0 indicates that PBLH remains level concerning the surface. When L = 1, the PBLH follows topography, in other words, higher terrains have higher PBLH values and vice-versa. Finally, L > 1 indicates that the amplitude (absolute value) of PBLH differences is larger than topographic ones, so higher (lower) terrain has much higher (lower) PBLH.
In the case of São Paulo city ∆z i = PBLH SEFAZ − PBLH IPEN and ∆z T = −57 m is the difference between the altitude of SEFAZ and IPEN lidar sites. To take advantage of the high spatial resolution of lidar systems (~7.5 m) used in São Paulo, cases L = 1 and L = 0 are replaced by L~1 and L~0, respectively. L~1 represents cases where L is closer to 1 than to 0, and PBLH tends to follow topography. As the difference between SEFAZ and IPEN altitudes ∆z T = −57 m, L~1 occurs when 57 m ≥ |∆z i | ≥ 29 m ( Table 1). On the other hand, L~0 indicates the cases where L is closer to 0 than to 1, and the PBLH tends to level when |∆z i | < 29 m. Table 1 summarizes how L and ∆z i are related in the case of São Paulo city.

Detection Algorithm of BB Events
Earlier studies show that in São Paulo city, the BB events are, in general, associated with the presence of high-concentration aerosol layers in the first 5 km of the atmosphere [10,25]. They occur between May and October when the number of forest fires in the center-west Brazil and Amazon regions and the burning of organic matter produced in the sugarcane crops in the countryside of São Paulo State [10,25] increases.
The intense aerosol loading combined with high solar radiation absorption capacity let these BB events be identified by the simultaneous increase in the daily values of Aerosol Optical Depth (AOD) and Aerosol Index (AI) [43,44]. Therefore, in this work, a new detection algorithm for BB events was developed and applied to São Paulo city. As displayed schematically in Figure 2, this algorithm combines key properties of lidar and satellite data. First, lidar data are used to find the aerosol layers. In the case of a positive outcome, daily values of AOD and AI are estimated for the lidar site and compared with corresponding daily values in the previous day.

PBLH Horizontal Variation
In this section, the diurnal evolution of differences between PBLH retrieved at IPEN (PBLHIPEN) and SEFAZ (PBLHSEFAZ) from lidar measurements is analyzed. Firstly, two case As there is no evidence that São Paulo city is affected by dust outbreaks, identifying the BB event is necessary only to verify whether daily values of AI remain positive during two consecutive days. To track the origin of the BB event is necessary to verify whether AOD daily values are increasing during these two consecutive days when AI values were positive. If both conditions are satisfied, it is plausible to infer that the aerosol layer detected in São Paulo was produced by a BB event and VIIRS data can be used to identify its origin.

PBLH Horizontal Variation
In this section, the diurnal evolution of differences between PBLH retrieved at IPEN (PBLH IPEN ) and SEFAZ (PBLH SEFAZ ) from lidar measurements is analyzed. Firstly, two case studies illustrate how the absence (23 July 2019) and presence (8 August 2019) of low clouds affect the relationship between PBLH and topographic difference during the day. Then, a statistical analysis of the PBLH behavior is discussed, considering observations carried out for 12 days selected is also presented. It is important to emphasize that sea-breeze circulation is observed with frequency in São Paulo city [31]; however, sea breeze circulation was not observed on 23 July and 8 August 2019.

Case 1: Absence of Low Clouds
The diurnal evolution of PBLH is shown (by black stars) in the RCS-Intensity curtain plots of 23 July 2019 at IPEN ( Figure 3a) and SEFAZ ( Figure 3b). Due to technical differences in the lidar systems and possible spatial variation in the aerosol concentrations, the RCS intensities observed in both sites are not the same. However, the vertical structure of the aerosol layers is quite similar. This day is characterized by the absence of low clouds and the aerosol layer adjacent to the surface is well mixed, consequently a well-defined PBL in both sites. There is a thin aerosol layer (purple dashed boxes in Figure 3a   At IPEN (Figure 3a), the PBLH remains practically constant from 09 to 11 LT, and it has a fast-growth stage from 12 to 14 LT reaching 1400.0 ± 7.5 m. Then, from 15 to 18 LT, PBLH remains practically constant and equal to 1452.0 ± 7.5 m.
At SEFAZ (Figure 3b), the PBLH displayed a longer fast-growth stage from 09 to 13 LT, reaching 1390.0 ± 7.5 m. At 15 LT, the PBLH reaches its maximum value of 1460.0 ± 7.5 m and maintains constant until 18 LT. Figure 4 display the diurnal evolution of PBLH SEFAZ -PBLH IPEN (∆z i ) and PBLH Levelness (L) on 23 July 2019. The higher differences occur during the growth stage because PBLH SEFAZ and PBLH IPEN displayed different growth rates (156.7 and 153.0 m·h −1 , respectively), which is probably associated with the land use of each region (IPEN-suburban, SEFAZ-urban). Comparatively, on most mornings, PBLH SEFAZ is higher and has larger growth rate values than PBLH IPEN . Such behavior is expected because the sensible heat flux in the urban region (SEFAZ) is expected to be higher than in the suburban region (IPEN). The maximum difference (∆z i = 300.0 ± 10.6 m) was observed at 12 LT when PBLH SEFAZ reached the mature stage and PBLH IPEN just started the fast-growth stage. After the fast growth stage, at 14 LT, PBLH IPEN and PBLH SEFAZ displayed similar behavior, reaching a remarkably similar maximum value and remaining almost constant until the end of the day. During this period, an intense reduction in the difference between the PBLH values occurs (∆z i = 8.0 ± 10.6 m).  At IPEN, PBLHIPEN displays a fast-growth stage from 09 to 15 LT, reaching its maximum value of 2600 ± 7.5 m. The presence of scattered low clouds affected the performance of the WTC method, which tends to retrieve PBLH as the base of the cloud [15]. During the afternoon, the presence of clouds became more frequent in both sites, and PBLH was retrieved below the actual PBLH, indicating that the cloud base was below the actual PBLH in the IPEN site. Between 16 and 18 LT, PBLH was retrieved below the actual PBLH, indicating that the cloud base was below the actual PBLH.
At SEFAZ, PBLHSEFAZ displays a fast-growth stage from 09 to 14 LT. Between 15:30 and 17 LT, the presence of scattered low clouds affected the performance of the PBLHdetection algorithm so that hourly values of PBLH retrieved from the WCT method were misplaced below the actual PBLH. At 18 LT, PBLHSEFAZ reaches its maximum daily value of 2400 ± 7.5 m. Regarding L, during the first two hours (09 to 10 LT), while the fast PBL-growth stage has not started, the PBLH tends to level (L~0). During the fast-growth stage (11 to 13 LT), PBLH varies horizontally opposite to topography because the growth rate at SEFAZ is higher than at IPEN (L < 0). However, at 14 LT, when both PBL are almost fully developed, the PBLH tends to level again (L~0). This result agrees with Stull [42], which indicates that the PBLH tends to level when it is fully developed.

Case 2: Presence of Low Clouds
The diurnal evolution of PBLH is shown (by black stars) in the RCS-Intensity curtain plots on 8 August 2019 at IPEN (Figure 5a) and SEFAZ (Figure 5b). This day is characterized by the presence of scattered low clouds at both sites.
At IPEN, PBLH IPEN displays a fast-growth stage from 09 to 15 LT, reaching its maximum value of 2600 ± 7.5 m. The presence of scattered low clouds affected the performance of the WTC method, which tends to retrieve PBLH as the base of the cloud [15]. During the afternoon, the presence of clouds became more frequent in both sites, and PBLH was retrieved below the actual PBLH, indicating that the cloud base was below the actual PBLH in the IPEN site. Between 16 and 18 LT, PBLH was retrieved below the actual PBLH, indicating that the cloud base was below the actual PBLH.
At SEFAZ, PBLHSEFAZ displays a fast-growth stage from 09 to 14 LT. Between 15 and 17 LT, the presence of scattered low clouds affected the performance of the PBL detection algorithm so that hourly values of PBLH retrieved from the WCT method w misplaced below the actual PBLH. At 18 LT, PBLHSEFAZ reaches its maximum daily va of 2400 ± 7.5 m.  At SEFAZ, PBLH SEFAZ displays a fast-growth stage from 09 to 14 LT. Between 15:30 and 17 LT, the presence of scattered low clouds affected the performance of the PBLHdetection algorithm so that hourly values of PBLH retrieved from the WCT method were misplaced below the actual PBLH. At 18 LT, PBLH SEFAZ reaches its maximum daily value of 2400 ± 7.5 m.
Like the previous case, PBLH IPEN and PBLH SEFAZ display different growth rates (266.7 and 166.7 m·h −1 , respectively). However, in this case, the PBLH IPEN grows faster than PBLH SEFAZ from 9 to 14 LT ( Figure 5), indicating that other causes rather than land use are affecting the diurnal evolution of PBLH in both sites.
Regarding L, it remains >1 during the fast-growth stage (09-13 LT) and most of the mature stage (15)(16)(17), indicating that the amplitude of the absolute values of PBLH differences are systematically larger than topographic ones. The higher ∆z i value (−500 ± 10.6 m) is observed during the mature stage period. Comparatively to the fastgrowth stage, the frequency of low clouds is higher during the mature stage in both sites ( Figure 5). At 14 LT, no clouds were identified by both lidars sites (Figure 5), the growth rate of PBLH IPEN and PBLH SEFAZ are equally attenuated, and the PBLH difference decreases significantly (∆z i = −48 ± 10.6 m). As indicated by L~1 at 14 LT (Figure 6), the PBLH tends to follow the topography. At 18 LT, the presence of low clouds was observed only in the IPEN region. As a result, PBLH differences become smaller than ∆z i (56 ± 10.6 m), and, as indicated by L < 0, PBLH difference varies in the opposite way to the topography one.
Based on the above analysis, it seems plausible to infer that the presence of clouds has affected the diurnal evolution of PBLH in both sites by keeping the PBLH differences between SEFAZ and IPEN negative, which contradicts the expected land-use effect ng and by increasing it in absolute terms during the afternoon.
stricting the vertical evolution of PBL to 1450 m in IPEN (Figure 3a) and SEFAZ ( Figure  3b). On 8 August 2019, the presence of low clouds disrupted the process of stabilization induced by the aerosol layers from BB events, allowing a much deeper vertical evolution of the PBL, but as indicated in the Figure 5a,b, clouds were more frequent over IPEN than SEFAZ, so that the PBL differences remain negative during almost the entire daytime and becomes more intense during the afternoon when the presence of low clouds are more intense at IPEN site.   (Figure 5b), higher aerosol layers remain above the PBL during the 17 LT. As will be shown in Section 4.2 below, these aerosol layers are composed of black carbons produced by BB events advected to São Paulo city during winter. During the daytime, they absorb incoming solar radiation, heating the atmosphere locally between 2 and 3 km. This localizing heat generates elevated inversion layers that, in turn, inhibit the vertical evolution of the PBL. This mechanism seems to be acting on 23 July 2019, when no low clouds were detected, restricting the vertical evolution of PBL to 1450 m in IPEN ( Figure 3a) and SEFAZ (Figure 3b). On 8 August 2019, the presence of low clouds disrupted the process of stabilization induced by the aerosol layers from BB events, allowing a much deeper vertical evolution of the PBL, but as indicated in the Figure 5a,b, clouds were more frequent over IPEN than SEFAZ, so that the PBL differences remain negative during almost the entire daytime and becomes more intense during the afternoon when the presence of low clouds are more intense at IPEN site.

All Campaign
The diurnal evolution of statistical properties of PBLH IPEN and PBLH SEFAZ , as well as their respective mean values (PBLH IPEN and PBLH SEFAZ ), retrieved from lidar performed in SEFAZ and IPEN during the 12 days of the field campaign, are presented in Figure 7. At 18 LT, PBLHSEFAZ has greater variability than PBLHIPEN and a predominance of smaller values. Such an effect is caused by the sea breeze, which systematically penetrates more than 50% of all days of the year in São Paulo city [31]. The sea breeze brings colder and more moist air from the coast, changing the vertical aerosol distribution and causing a high aerosol concentration in the lower part of the PBL. Then, methods based on the vertical aerosol gradient (e.g., WCT) tend to underestimate the PBLH in this situation [7]. Analyzing high-resolution data from numerical simulation with WRF Model, Ribeiro et al. [32] also show that the passage of sea breeze during the afternoon induces a significant drop in the PBLH in São Paulo city.  . From 09 to 10 LT, the average zi remains below |28| m so that during this period, the PBLH tends to level (L ~ 0). At 11 LT, PBLHSEFAZ display a large growth rate, zi intensifies, reaching 55 ± 30.7 m and the average horizontal variation of PBLH is opposite to the topography (L < 1). Between 12 and 13 LT, the average At 18 LT, PBLH SEFAZ has greater variability than PBLH IPEN and a predominance of smaller values. Such an effect is caused by the sea breeze, which systematically penetrates more than 50% of all days of the year in São Paulo city [31]. The sea breeze brings colder and more moist air from the coast, changing the vertical aerosol distribution and causing a high aerosol concentration in the lower part of the PBL. Then, methods based on the vertical aerosol gradient (e.g., WCT) tend to underestimate the PBLH in this situation [7]. Analyzing high-resolution data from numerical simulation with WRF Model, Ribeiro et al. [32] also show that the passage of sea breeze during the afternoon induces a significant drop in the PBLH in São Paulo city. Figure 8 presents the diurnal evolution of the levelness parameter and the difference between PBLH SEFAZ and PBLH IPEN (∆z i ). From 09 to 10 LT, the average ∆z i remains below |28| m so that during this period, the PBLH tends to level (L~0). At 11 LT, PBLH SEFAZ display a large growth rate, ∆z i intensifies, reaching 55 ± 30.7 m and the average horizontal variation of PBLH is opposite to the topography (L < 1). Between 12 and 13 LT, the average PBLH displays the most intense growth rate in both sites and ∆z i becomes greater, so that the PBLH varies horizontally, amplifying topographic differences (L > 1). From 14 to 16 LT, the period in which PBLH SEFAZ and PBLH IPEN have little variation, ∆z i is again smaller than 28 m, and PBLH tends to levelling (L~0). At 17 LT, due to the reduction in the mean value of PBLH IPEN , PBLH varies horizontally as opposed to topography (L < 0). Then, at 18 LT, the PBLH SEFAZ reduces sharply so that the PBLH varies, intensifying the topographical differences (L > 1). Such a phenomenon occurs due to the presence of the sea breeze, which arrives first at the SEFAZ, later spreading to the IPEN region. However, due to slow sea breeze displacement in urban areas [32,42], at 18 LT, only the SEFAZ region was reached by it. So, PBLH SEFAZ has a significant drop (Figure 7), becoming lower than PBLH IPEN .
Atmosphere 2022, 13, x FOR PEER REVIEW 11 of 17 due to slow sea breeze displacement in urban areas [32,42], at 18 LT, only the SEFAZ region was reached by it. So, PBLHSEFAZ has a significant drop (Figure 7), becoming lower than PBLHIPEN. It is important to emphasize that the effects of thermal stabilization associated with the presence of aerosol layers associated with BB events on average seem to be canceled out by the effect of low clouds so that PBLH differences between SEFAZ and IPEN are very small. Figure 9 shows the RCS-Intensity curtain plots measured by the lidar at the SEFAZ It is important to emphasize that the effects of thermal stabilization associated with the presence of aerosol layers associated with BB events on average seem to be canceled out by the effect of low clouds so that PBLH differences between SEFAZ and IPEN are very small.  Figure 9 shows the RCS-Intensity curtain plots measured by the lidar at the SEFAZ on 19 August 2019. From 06 to 13 LT, an intense aerosol plume (dashed violet box) extends from 1000 m to the PBLH SEFAZ . In the presence of intense aerosol plumes, such as the BB event on 19 August, the lidar methods based on the aerosol gradient tend to retrieve PBLH at the top of the aerosol layer. Therefore, the behavior of PBLH SEFAZ between 09 and 12 LT does not necessarily represent the real top of the PBL [2,13]. Indeed, on average, the PBLH displays a fast-growth stage during this period (Figure 7b). From 22 LT until the end of the measurements at 24 LT, a shallow layer with high aerosol concentration appeared near the surface, causing an intense backscattering and impeding laser signal to reach above this shallow layer. It is important to emphasize that the effects of thermal stabilization associated with the presence of aerosol layers associated with BB events on average seem to be canceled out by the effect of low clouds so that PBLH differences between SEFAZ and IPEN are very small. Figure 9 shows the RCS-Intensity curtain plots measured by the lidar at the SEFAZ on 19 August 2019. From 06 to 13 LT, an intense aerosol plume (dashed violet box) extends from 1000 m to the PBLHSEFAZ. In the presence of intense aerosol plumes, such as the BB event on 19 August, the lidar methods based on the aerosol gradient tend to retrieve PBLH at the top of the aerosol layer. Therefore, the behavior of PBLHSEFAZ between 09 and 12 LT does not necessarily represent the real top of the PBL [2,13]. Indeed, on average, the PBLH displays a fast-growth stage during this period (Figure 7b). From 22 LT until the end of the measurements at 24 LT, a shallow layer with high aerosol concentration appeared near the surface, causing an intense backscattering and impeding laser signal to reach above this shallow layer.          Inversion data from SP-EACH AERONET sunphotometer station were employed to derive the log-normal size distribution of aerosol from 17 to 19 August 2019. As shown in Figure 12, the log-normal distribution on 17 August changes during the day. However, Inversion data from SP-EACH AERONET sunphotometer station were employed to derive the log-normal size distribution of aerosol from 17 to 19 August 2019. As shown in Figure 12, the log-normal distribution on 17 August changes during the day. However, fine mode aerosols are dominant during the whole period and there is strong evidence that most aerosols come from the BB event produced by the wildfires on the Brazil-Bolivia frontier.  ( ) ( ) ⁄ is the particle volume density and is the particle radius.

Biomass Burning Detection
According to Cazorla et al., 2012 [45], the chemical composition of aerosol can be identified from the Ångström Matrix. It consists of a dispersion diagram of Scattering Ångström Exponent (SAE) by Absorption Ångström Exponent (AAE). The sub-sections of this diagram define aerosol composition by grouping absorbing aerosol types [46,47]. During the aerosol transport event that occurred between 17 and 19 August over MRSP, the Ångström Matrix in Figure 13 shows the presence of small particles with low absorption, with some cases of black carbon aerosol composition. . dV(r)/dln(r) is the particle volume density and r is the particle radius.
According to Cazorla et al., 2012 [45], the chemical composition of aerosol can be identified from the Ångström Matrix. It consists of a dispersion diagram of Scattering Ångström Exponent (SAE) by Absorption Ångström Exponent (AAE). The sub-sections of this diagram define aerosol composition by grouping absorbing aerosol types [46,47]. During the aerosol transport event that occurred between 17 and 19 August over MRSP, the Ångström Matrix in Figure 13 shows the presence of small particles with low absorption, with some cases of black carbon aerosol composition. Figure 12. Aerosol size distribution was retrieved from AERONET sunphotometer measurements at the EACH station in São Paulo city (23°48′ S, 46°49′ W, 754 m asl) on 17 August 2019. Solid circlelines numbers in the legend (top-left) indicate local time, particle radius (µm), number, and volume density (µm 3 µm −2 ).
( ) ( ) ⁄ is the particle volume density and is the particle radius.
According to Cazorla et al., 2012 [45], the chemical composition of aerosol can be identified from the Ångström Matrix. It consists of a dispersion diagram of Scattering Ångström Exponent (SAE) by Absorption Ångström Exponent (AAE). The sub-sections of this diagram define aerosol composition by grouping absorbing aerosol types [46,47]. During the aerosol transport event that occurred between 17 and 19 August over MRSP, the Ångström Matrix in Figure 13 shows the presence of small particles with low absorption, with some cases of black carbon aerosol composition.  Therefore, based on the methodology presented in Section 3.3, it is possible to conclude that the intense aerosol layer observed on 19 August 2019 represents a BB event. Such a result is in accordance with Pereira et al. (2021) [48], who from a combination of back trajectory analyses (from Hybrid Single-Particle Lagrangian Integrated Trajectory [49]) and images from satellite (GOES-16 Advanced Baseline Instrument and Moderate Resolution Imaging Spectroradiometer) observed air masses transporting plumes of biomass burning aerosols from Bolivia and the Amazon Basin to São Paulo state. Such a BB event was widely publicized in the press due to its unprecedented effects, such as "black rain" [48], which could be detected in several regions of the São Paulo state.
In the 12 days of measurements, strong aerosol plumes above or on the PBL top were found in 3 of them (23 July, 10 and 19 August 2019), but only on 10 August variations in AOD and AI could be associated with the presence of material from biomass burning.

Conclusions
Elastic lidar systems have been widely applied in studies related to PBL structure. The high spatial and temporal resolutions of this type of remote sensing system allow it to be used in several environmental applications. In this work, data from two elastic LIDAR systems, which operated simultaneously at approximately 10.7 km of distance, were used to analyze the homogeneity of the PBLH in the city of São Paulo. In addition, a new procedure based on elastic lidar and satellites data were used to identify the impact of BB events in the aerosol layer detected above the PBL.
Except for the cases with the presence of sea breezes, in general, higher differences in the PBLH horizontal homogeneity (around 100 m) are observed during its growth period. Before the PBL growth stage and when the PBLH is fully developed, it was observed that PBLH tends to levelling (L~0). Sea breeze moves slowly in urban areas, so its presence makes the PBLH vary horizontally, amplifying topographical differences (L > 1). Furthermore, a stabilization of the atmosphere caused by the presence of aerosol layers associated with BB events was observed. However, such a result needs a more intensive investigation using radiosonde data to identify the presence of thermal inversion induced by these aerosol layers.
In addition, the synergy between the elastic lidar and satellites data enabled the detection of a BB event, as well as identifying its origin. Such a result demonstrates that the proposed methodology is an efficient and easy-applicable tool to detect BB events in São Paulo city. However, it is necessary to extend the application to a larger set of cases, using AE as a benchmark.
Therefore, these results reinforce the great importance and applicability of elastic lidar systems in PBL studies. So that the combination of a set of lidar systems or a synergy use of lidars and other remote sensing systems can provide a better understanding of certain phenomena.