Anthropogenic Aerosol Contribution to Urban Environments. Aerosol Optical Depth Analysis During COVID Pandemic and Post-Pandemic Periods

1. Introduction

Since the Great London Fog, the influence of meteorological conditions, air quality, and pollutant emissions on public health has received sustained attention. Despite advances in air-quality management, rapid urbanization and industrialization continue to make urban air pollution a significant and persistent challenge (Wu et al., 2021; Via et al., 2021; Eleftheriadis et al., 2021; Mananga et al., 2023). Of particular concern is fine particulate matter with an aerodynamic diameter of ≤2.5 μm (PM2.5). These particles can penetrate deep into the respiratory system and enter the bloodstream, increasing the risk of respiratory and cardiovascular diseases, particularly among vulnerable populations (Mananga et al., 2023). In New York City, PM2.5 exposure from all sources is estimated to contribute annually to approximately 2,000 deaths and 5,150 emergency department visits and hospitalizations related to respiratory and cardiovascular disease (NYC Environmental & Health Data Portal, 2021). Similarly, a study in Rome reported that a 10 μg m⁻³ increase in PM2.5 concentration was associated with a 4% increase in the risk of nonaccidental mortality (Cesaroni et al., 2013).

This study evaluates aerosol and air-quality conditions over New York City during June–July 2020, representing the acute phase of the COVID-19 pandemic, and the corresponding period in 2023, representing post-pandemic conditions. The intervening years were excluded because partial pandemic restrictions remained in place. The effects of pandemic-related changes in anthropogenic activity on particulate pollution remain debated. Whereas Zangari et al. (2020) reported no significant changes in particulate pollution, Sarfraz et al. (2020) and Shehzad et al. (2021) identified substantial variations. These contrasting findings demonstrate the need for rigorous, localized analyses to quantify the influence of anthropogenic activity on urban air quality.

To characterize aerosol conditions and transport processes, this study employs a suite of complementary observations, including a Lufft CHM15k ceilometer. Lidar remote-sensing techniques provide high-resolution vertical information on aerosol distributions and transport pathways and are widely used for atmospheric particulate-matter monitoring (Vladutescu et al., 2012). Spaceborne lidar observations from CALIPSO, for example, have documented the horizontal and vertical evolution of pollution plumes over major urban areas, including Mexico City (Kar et al., 2015), while ground-based Doppler lidars have resolved fine-scale turbulence and pollutant dispersion in densely populated coastal environments (Arend et al., 2020). Multiwavelength elastic–Raman lidar techniques have further demonstrated the ability to retrieve aerosol and cloud properties simultaneously (Wu et al., 2011). Here, ceilometer-derived aerosol observations within the planetary boundary layer are used to assess the influence of human activity on air quality in New York City.

Lidar retrievals were performed using the Fernald algorithm, a widely applied inversion method for deriving aerosol extinction profiles from elastic-backscatter signals. The method solves the lidar equation by assuming a constant or parameterized lidar ratio and applying a boundary condition at a selected reference range, thereby estimating the vertical distribution of aerosols (Fernald et al., 1984). Building on the earlier work of Klett (1981), the Fernald approach provides a practical framework for aerosol retrieval under varying atmospheric conditions.

Interpretation of lidar-derived aerosol profiles requires consideration of planetary boundary layer (PBL) dynamics and chemistry. The planetary boundary layer height (PBLH) regulates the vertical distribution, mixing, and dispersion of air pollutants in urban environments. Relative to rural areas, urban PBL evolution is influenced by the urban heat-island effect, enhanced surface roughness, and anthropogenic heat emissions (Kim et al., 2021). These factors can promote a deeper and more turbulent daytime boundary layer, enhancing vertical mixing and temporarily diluting near-surface pollutant concentrations. In contrast, nighttime stabilization often produces a shallow PBL that confines pollutants near the surface and increases exposure risk (Jacobson, 2002). The PBLH exhibits pronounced diurnal variability and responds to wind speed, solar radiation, and cloud cover. Consequently, characterization of PBL dynamics is essential for interpreting aerosol lidar profiles and assessing urban air quality. Recent work has further demonstrated that the vertical distribution of wildfire smoke is strongly coupled to PBL behavior, including nocturnal boundary-layer collapse and daytime convective mixing (Huang et al., 2024).

Local emissions from road traffic, shipping, and aviation, as well as aerosols transported over long distances, can influence urban aerosol loading. Under stagnant conditions, locally emitted pollutants may dominate near-surface concentrations, whereas regional transport can introduce dust, wildfire smoke, and other aerosol types from distant source regions (Vladutescu et al., 2013). Atmospheric stratification can maintain elevated aerosol layers above the PBL; subsequent boundary-layer growth and entrainment may affect surface-level pollutant concentrations. Fine particles, including PM2.5 and black carbon, may therefore contribute to air-quality degradation through both local emissions and long-range transport processes (Hastie et al., 1993; Schumann, 1996).

The COVID-19 pandemic created an unprecedented opportunity to examine aerosol optical properties and atmospheric pollution under substantially reduced anthropogenic activity. Previous studies have reported that quarantine measures reduced emissions and improved urban air quality (Eleftheriadis et al., 2021). Consistent with these observations, the present study identifies aerosol-related differences of approximately a factor of two to three between the pandemic and post-pandemic periods. Although New York City’s initial stay-at-home order began on 20 March 2020, sparse data availability during March–May motivated the selection of June–July 2020 for the pandemic analysis and the corresponding months of 2023 for the post-pandemic comparison. New York City entered the first phase of reopening on 8 June 2020, following fulfillment of the State’s reopening criteria; occupancy restrictions and other safety protocols remained in effect during the study period.

The remainder of this paper is organized as follows. Section 2 describes the instrumentation and data sets. Section 3 presents the methodology and theoretical background. Section 4 provides the results and discussion, and Section 5 summarizes the main conclusions. The data used in this study are available from the City College of New York Data Server, NASA AERONET, the New York State Department of Environmental Conservation air-monitoring network, and the University of Wyoming atmospheric-sounding archive.

2. Instruments

The analysis integrates observations from a Lufft CHM15k ceilometer lidar, a CIMEL sun photometer, radiosondes, a 1405-DF TEOM ambient PM2.5 monitoring system, and HYSPLIT backward-trajectory simulations. These complementary data sets were used to characterize aerosol optical properties, vertical aerosol structure, local particulate-matter concentrations, and potential aerosol transport pathways.

2.1. Lufft CHM15k Ceilometer

The Lufft CHM15k is an autonomous, single-wavelength elastic-backscatter ceilometer lidar operating at 1064 nm (Wiegner et al., 2012, 2019). It measures range-resolved atmospheric backscatter signals to heights of up to 15 km, with a vertical resolution of 5–15 m, depending on the selected measurement range. The temporal resolution is configurable from 15 to 60 s. In addition to raw backscatter profiles, the instrument provides cloud-base height for up to nine cloud layers, penetration depth, maximum detection range, vertical visibility, sky-condition classification, and cloud amount.

The ceilometer comprises a laser transmitter, receiver, optical system, environmental sensors, and data-acquisition unit. The transmitter emits laser pulses with durations of 1–5 ns at a pulse repetition rate of 5–7 Hz and pulse energies of 7–9 μJ. The emitted beam is directed into the atmosphere through the optical system, and the receiver records the intensity of the radiation backscattered by atmospheric constituents as a function of time, from which altitude-resolved backscatter profiles are derived. Internal and external temperature sensors monitor instrument operating conditions. In this study, CHM15k backscatter observations were used to investigate cloud and aerosol vertical distributions, identify elevated plumes, and estimate the evolution of the planetary boundary layer.

2.2. CIMEL CE318-T Sun Photometer

The CIMEL CE318-T sun photometer measures direct solar irradiance at multiple wavelengths using spectral filters and silicon and InGaAs detectors. These measurements are used to derive aerosol optical depth (AOD), the Ångström exponent, and precipitable water vapor (PWV). The instrument deployed at the City College of New York (CCNY) campus is part of the NASA Aerosol Robotic Network (AERONET) (Holben et al., 1998) and is shown in Figure 1b.

The CE318-T consists of a control unit, optical sensors, a collimator, a robotic sun-tracking system, and a tripod. The collimator and robotic tracker maintain accurate alignment with the Sun to enable direct-sun irradiance measurements. The tracker operates over azimuth and zenith ranges of 0–360° and 0–180°, respectively, with an angular resolution of 0.003°. The instrument is equipped with nine spectral channels centered at 340, 380, 440, 500, 675, 870, 937, 1020, and 1640 nm, spanning the ultraviolet, visible, and near-infrared spectral regions. Solar irradiance measured sequentially through these filters is transmitted to the data-acquisition and processing unit, from which column-integrated aerosol and water-vapor properties are retrieved.

2.3. Radiosondes

Radiosonde observations were used to characterize the vertical profiles of atmospheric temperature, relative humidity, pressure, wind speed, and wind direction. Radiosondes are carried aloft by weather balloons, also known as sounding balloons, from the surface to the upper atmosphere. A typical radiosonde includes temperature, humidity, and pressure sensors, a Global Positioning System (GPS) receiver, and a transmitter. Operational upper-air soundings are commonly launched twice daily, at 00:00 and 12:00 UTC, providing globally consistent measurements of atmospheric structure. In this study, radiosonde data were obtained from the University of Wyoming Upper Air Soundings archive. The platform provides historical and near-real-time sounding data for stations worldwide and enables the selection of data by geographic region, station, date, and launch time. A representative radiosonde system is shown in Figure 1d.

2.4. TEOM 1405-DF Ambient Particulate-Matter Monitor

Ambient particulate-matter concentrations were measured using a Thermo Scientific Tapered Element Oscillating Microbalance (TEOM) 1405-DF system (Figure 1c). The TEOM 1405-DF provides continuous, near-real-time measurements of particulate-matter mass concentration. Depending on the sampling configuration, its dual-filter design can measure PM2.5 and either coarse particulate matter (PM10–2.5) or other particle-size fractions.

The instrument draws ambient air through a filter mounted on a tapered oscillating element. Particle accumulation on the filter increases the element mass and produces a measurable decrease in oscillation frequency. The measured frequency shift is converted to particulate-matter mass concentration, enabling high-temporal-resolution monitoring. Temperature and humidity control systems reduce measurement artifacts associated with water vapor and volatile particulate components. The TEOM 1405-DF is widely used for regulatory air-quality monitoring, long-term pollution-trend assessment, and exposure-related studies.

3. Methodology

The ceilometer measures the backscatter intensity of light in the time domain to determine cloud base height and aerosol backscatter profile. In this retrieval, we assumed a suitable LIDAR ratio and used the radiosonde data. Based on these, we find the extinction and backscatter coefficient by using Fernald’s method summarized below (Fernald, F. G., 1984). Further, we use the extinction coefficient to calculate the aerosol optical depth. Finally, we compare the difference in AOD between our results and the AERONET data. Hence, we apply the above procedures for two summer months to compare the differences between 2020 and 2023. Furthermore, we use the HYSPLIT MODEL for trajectory analysis, and CIMEL sunphotometer derived AOD to cross-validate and better understand the ceilometer’s measurements. This multi-instrument approach helps us verify and improve our understanding of atmospheric composition.

3.1. LIDAR Ratio, Backscatter, and Extinction

The LIDAR ratio is a crucial parameter in remote sensing and atmosphere studies. The LIDAR ratio is defined as the ratio of the particle extinction coefficient (α) and the particle backscattering coefficient (β) (Weitkamp, 2005). In our research, we assumed a LIDAR ratio of 40 sr at 1064-nm , a commonly used value for anthropogenic aerosols (Li, D. 2022; Omar et al., 2009).

The extinction coefficient is the sum of loss from absorption and light scattering. It refers to the reduction of the intensity of light as it propagates through the atmosphere due to scattering and absorption by particulates and gases. It quantifies how aerosols and gases affect the transmission of light, leading to phenomena such as haziness, reduced visibility, and color changes in the sky. High extinction values can indicate the presence of dense aerosols, clouds or pollutants. This presence effectively attenuates the LIDAR signal by absorption or scattering. In Eq. (1), the extinction coefficient is the sum of the scattering by molecules, scattering by particles, absorption by molecules, and absorption by particles.

$$\begin{array}{c}{\alpha }_{ \lambda }\left(Z\right)={\alpha }_{s,mol , \lambda }\left(Z\right)+{\alpha }_{s, par, \lambda }\left(Z\right)+{\alpha }_{a, mol, \lambda }\left(Z\right)+{\alpha }_{a, par, \lambda }\left(Z\right) \end{array}$$

(1)

Backscatter refers to the light that is reflected directly back to the LIDAR sensor after it interacts with atmospheric particles or molecules. It provides information about the properties and concentrations of aerosols and clouds. The strength and characteristics of the backscattered light help in determining particulate distribution and concentration, as well as in evaluating atmospheric conditions such as visibility and pollution levels. High backscatter values suggest the presence of larger or denser concentrations of particles, whereas lower backscatter values may indicate cleaner air or the presence of smaller particles. In Eq. (2), the backscatter coefficient is the sum of the backscatter by molecules and particles.

$$\begin{array}{c}{\beta }_{ \lambda }\left(Z\right)={\beta }_{mol , \lambda }\left(Z\right)+{\beta }_{par, \lambda }\left(Z\right) \end{array}$$

(2)

The Fernald method is a commonly used technique in LIDAR remote sensing to derive profiles of atmospheric backscatter and extinction coefficients from measured LIDAR signals. It assumes a known relationship between the backscatter and extinction (usually the LIDAR ratio) and requires a reference point in the atmosphere where the backscatter coefficient is known. Fernald explained the LIDAR equation, P(Z) is the return signal that is proportional to the received power from a scattering volume at a slant range Z as shown in Eq. (3)

$$\begin{array}{c}P\left(Z\right)=EC{Z}^{-2}\left[{\beta }_{par}\left(Z\right)+{\beta }_{mol}\left(Z\right)\right]T_{par}^2\left(Z\right)T_{mol}^2\left(Z\right) \end{array}$$

(3)

Where E is the output energy pulse that is proportional to the transmitted energy; C is the calibration constant of the instrument, which includes losses in the transmitting and receiving optics and the effective receiver aperture; β_par(Z) and β_mol(Z) are the backscattering cross sections of the aerosols and molecules at slant range Z, respectively. In this equation, $T_{par}\left(Z\right)=exp\exp \left[-{\int }_0^Z{\alpha }_{par}\left(Z\right)dZ\right]$, is the aerosol transmittance, and $T_{mol}\left(Z\right)=exp\exp \left[-{\int }_0^Z{\alpha }_{mol}\left(Z\right)dZ\right]$, is the molecular atmosphere transmittance. Additionally, ${\alpha }_{par}$ (Z) and ${\alpha }_{mol}$ (Z) are the extinction coefficients of the aerosols and molecules at range Z, respectively. (Fernald et al. 1984)

According to the Fernald inversion method (Fernald, 1984), the aerosol backscatter coefficient is retrieved recursively from a selected reference altitude by solving the lidar equation in discrete form. Let $X\left(Z\right)=P\left(Z\right)Z^2$ denote the range-corrected signal and △Z the range resolution. A reference altitude Z(I) is selected in a relatively clean atmospheric region, where the aerosol backscatter coefficient is specified using an assumed scattering ratio as the boundary condition (here I represents the current iteration layer). Eq. (4) incorporates the effects of aerosol extinction between adjacent data points, range ΔZ apart. In this case is shown for the upper range, with the lower range being calculated similarly to determine $A(I-1,I$). Starting from the reference point $Z\left(I\right)$, the total backscatter coefficient at the adjacent lower range gate $Z\left(I-1\right)$is calculated using Eq. (5), while the value at the adjacent upper range gate $Z\left(I+1\right)$is calculated using Eq. (6). In these equations, ${\beta }_{par}$and ${\beta }_{mol}$represent the aerosol and molecular backscatter coefficients, respectively, and the exponential terms account for the molecular transmission correction between adjacent range gates.

After the aerosol backscatter coefficient is obtained, the corresponding aerosol extinction coefficient is derived using the assumed aerosol lidar ratio. The extinction-form recursive equations are given in Eq. (7) and (8), where $S_{par}$and $S_{mol }$are the aerosol and molecular lidar ratios, taken as 40 sr and 8*3.14/3 sr, respectively. Eq. (7) gives the extinction-related solution for the lower range gate $Z\left(I-1\right)$, while Eq. (8) gives the corresponding solution for the upper range gate $Z\left(I+1\right)$. This backward and forward recursive procedure is applied to each individual time profile to retrieve the time-height distribution of aerosol backscatter and extinction coefficients.

$$\begin{array}{c}A\left(I,I+1\right)=\left(S_{par}-S_{mol}\right)\left[{\beta }_{mol}\left(I\right)+{\beta }_{mol}\left(I+1\right)\right]\Delta Z\end{array}$$

(4)

$$\begin{array}{c}{\beta }_{par}\left(I-1\right)+{\beta }_{mol}\left(I-1\right)={\displaystyle \frac{X\left(I-1\right)exp \left[+A\left(I-1, I\right)\right] }{\frac{X\left(I\right)}{{\beta }_{par}\left(I\right)+{\beta }_{mol}\left(I\right)}+S_{par}\left[X\left(I\right)+X\left(I-1\right)exp \left[+A\left(I-1, I\right)\right] \right]△Z}}\end{array}$$

(5)

$$\begin{array}{c}{\beta }_{par}\left(I+1\right)+{\beta }_{mol}\left(I+1\right)={\displaystyle \frac{X\left(I+1\right)exp \left[-A\left(I, I+1\right)\right] }{\frac{X\left(I\right)}{{\beta }_{par}\left(I\right)+{\beta }_{mol}\left(I\right)}-S_{par}\left[X\left(I\right)+X\left(I+1\right)exp \left[-A\left(I, I+1\right)\right] \right]△Z}} \end{array}$$

(6)

$$\begin{array}{c}{\alpha }_{par}\left(I-1\right)+{\displaystyle \frac{S_{par}}{S_{mol}}}{\alpha }_{mol}\left(I-1\right)={\displaystyle \frac{X\left(I-1\right)expexp \left[+A\left(I-1, I\right)\right] }{\frac{X\left(I\right)}{{\alpha }_{par}\left(I\right)+\frac{S_{par}}{S_{mal}}{\alpha }_{mol}\left(I\right)}+\left[X\left(I\right)+X\left(I-1\right)expexp \left[+A\left(I-1, I\right)\right] \right]△Z}} \end{array}$$

(7)

$$\begin{array}{c}{\alpha }_{par}\left(I+1\right)+{\displaystyle \frac{S_{par}}{S_{mol}}}{\alpha }_{mol}\left(I+1\right)={\displaystyle \frac{X\left(I+1\right)expexp \left[-A\left(I, I+1\right)\right] }{\frac{X\left(I\right)}{{\alpha }_{par}\left(I\right)+\frac{S_{par}}{S_{mal}}{\alpha }_{mol}\left(I\right)}-\left[X\left(I\right)+X\left(I+1\right)expexp \left[-A\left(I, I+1\right)\right] \right]△Z}} \end{array}$$

(8)

To satisfy the inputs of Fernald method, the molecular backscatter is necessary (Geisinger et al. 2016). The backscatter coefficient by molecules is equal to the molecule number density (Eq. (9)) multiplied by the differential backscatter cross-section as in Eq. (10).

The molecule number density is equal to the pressure in Pascal divided by the product of Boltzmann constant k (k = 1.380649 × 10⁻²³ m² kg s⁻² K⁻¹) and temperature in Kelvin (see Eq. (9)). The pressure and the temperature were retrieved from the BNL (90 km east of New York City) radiosonde readings available on the Atmospheric Soundings website at weather.uwyo.edu (University of Wyoming Atmospheric Soundings website).

$$\begin{array}{c}{N}_{mol , \lambda }\left(Z\right)={\displaystyle \frac{P\left(Z\right)}{k*T\left(Z\right)}} \end{array}$$

(9)

This molecular number density was used in Eq. (10) to solve for the molecular backscatter coefficient ${\beta }_{mol , \lambda }\left(Z\right)$. Because only two radiosondes measurements per day were available (00 and 12 UTC), the atmospheric profiles were interpolated to provide continuous input for the retrieval calculations.

$$\begin{array}{c}{\beta }_{mol , \lambda }\left(Z\right)=N_{mol , \lambda }\left(Z\right)\left({\displaystyle \frac{d{\sigma }_{s, mol , \lambda }}{d\Omega }}\right) \end{array}$$

(10)

Where ${\sigma }_{s, mol , \lambda }$ is the scattering cross section of molecules at the given wavelength which according to Bucholtz’ look-up table (Bucholtz et al., 1995) is 2.734* 10⁻²⁸ cm² at wavelength 1100 nm (closest to the ceilometer wavelength of 1064 nm).

The Ceilometer calibration constant was implemented following the method proposed by O’Connor et al. (2004). In this procedure, an optically thick liquid cloud was selected as the calibration target. The liquid cloud lidar ratio (extinction/backscatter), $S=18 sr$, and the multiple-scattering correction factor, with values between 0.5 (thick optical clouds) and 1, in our case this factor is $\eta =0.65$. Under the optically thick cloud assumption, the theoretical integrated attenuated backscatter,${ \beta }_{theory}$, can be expressed as Eq. (11) and represents the expected integrated backscatter for an opaque liquid cloud.

$$\begin{array}{c}{\beta }_{theory}={\displaystyle \frac{1}{2\eta S}} \end{array}$$

(11)

$$\begin{array}{c}\beta ={\int }_{z_{base}}^{z_{top}}{\beta }_{meas}\left(z\right)\hspace{0.17em}dz \end{array}$$

(12)

The measured integrated backscatter, β, was obtained by vertically integrating the measured range-corrected backscatter signal between the detected cloud base and cloud top and is given in Eq. (12). The calibration coefficient was then calculated as the ratio between the measured integrated backscatter and the theoretical integrated backscatter, as indicated in Eq. (13).

$$\begin{array}{c}C={\displaystyle \frac{\beta }{{\beta }_{theory}}} \end{array}$$

(13)

3.2. Aerosol Optical Depth

Aerosol optical depth is a measure of the extinction of solar radiation by aerosol particles in the atmosphere. It is defined as the vertical integration of the extinction coefficient.

$$\tau =\int \alpha \left(Z\right)dZ$$

(14)

An integrated AOD over CCNY is available through the NASA AERONET website for the CIMEL (Vladutescu et al. 2013), and retrieved from extinction profiles measured by the LIDAR.

3.3. NOAA HYSPLIT Trajectory Model

The Hybrid Single-Particle LaGrangeian Integrated (Trajectory HYSPLIT) model tracks movement of air masses either forward or backward. We used longitude and latitude for our location, 40.82131° N and 73.94904° W. Then, we used HYSPLIT to track the aerosol origins. Some important factors are altitude, rainfall, mixing layer height, etc. Therefore, with the help of this model we can track the aerosol movement into and out of the zone of interest.

The results and analysis of the measurements with the above instruments are given in the section 4, below.

4. Results and Analysis

In this section, we present some of the LIDAR vertical profile used in the analysis of the summer months of June and July of 2020 and 2023. We also validated the AOD retrievals using the CIMEL sunphotometer, and compared the backscatter signal to the PM2.5 measurements, and looked into the possible sources.

4. a. LIDAR Retrievals

The vertical profiles of calibrated range-corrected return signals for the months of June-July 2020 are illustrated in Figure 2. These profiles indicate that the signal between 15–24 UTC hours are strongly contaminated by solar background above the PBL in clear-air regions, requiring careful daily analysis and selection of a reference level with low to no aerosol load for accurate retrievals. The colorbar values indicate the attenuated backscatter coefficient (/km/sr) or calibrated range-corrected return signals.

The recorded lidar signals were further processed to retrieve vertical profiles of attenuated backscatter and derived aerosol extinction coefficients, as shown in Figure 3 and Figure 4. These figures illustrate aerosol-loading conditions observed during the 2020 and 2023 measurement periods. For each year, representative cases spanning relatively high aerosol loading to comparatively clean atmospheric conditions were selected: 10 June, 22 June, and 27 July 2020, and 1 June, 6 July, and 22 July 2023. Instrument power calibrations were corrected for the respective measurement periods. These cases demonstrate typical temporal variations in aerosol and cloud conditions.

As described in Section 3, the retrieval is based on the Fernald iterative inversion method, which requires selection of an appropriate molecular reference altitude. An unsuitable reference altitude can introduce artifacts and noise into the retrieved extinction and backscatter profiles. Therefore, reference regions were selected daily under relatively stable atmospheric conditions. Cloud-contaminated profiles were excluded before the extinction retrieval because clouds strongly attenuate and scatter the lidar signal, rendering aerosol retrievals within and below cloud layers unreliable. Cloud detections are identified by purple markers in the backscatter-profile figures. Figure 3 presents cases of enhanced aerosol backscatter (Figure 3a), moderate aerosol loading (Figure 3b), and relatively low aerosol loading (Figure 3c). A common color scale is used in all panels to enable direct comparison.

Sample selected days after the pandemic are illustrated in Figure 4 for June 1^st, 2023 (Thursday), July 6^th, 2023 (Thursday), and July 22^nd, 2023 (Saturday).

The comparison of LIDAR retrieved AOD for 2023, indicates increased particulate pollution. We extended the comparison to cover two months (June and July). The LIDAR based aerosol optical depth was calculated for the entire duration of June and July of the two years and plotted for comparison in Figure 5. The retrieval approach is based on the methodology detailed in section three of the paper.

In this analysis, we applied a running average filter to remove noise, smoothen the strong transition between day and night, and to make it easier to observe trends in the data. The AOD relative difference was calculated based on the following formula:

$$\begin{array}{c}{AOD}_{re{l}_{diff}}={\displaystyle \frac{{AOD}_{2023}-{AOD}_{2020}}{{AOD}_{2020}}}\end{array}$$

(15)

Aerosol optical depth (AOD) values during summer 2023 were observed to be significantly higher than those observed during the corresponding period in 2020. Several elevated AOD episodes in June and July 2023 were associated with enhanced aerosol loading, including smoke transported from Canadian wildfire events. Some days identified as statistical outliers (see Section V) were also affected by intense precipitation, including 22 June 2023. Backscatter profiles for representative cases are presented in Figure 6. The profile for 7 June 2023 (Figure 6a) shows enhanced aerosol backscatter throughout much of the planetary boundary layer (PBL), consistent with the presence of transported wildfire smoke. In contrast, on 22 June, the ceilometer detected extensive cloud cover and rainfall; magenta markers indicate cloud-base altitude as determined by the proprietary Lufft software. The corresponding AOD time series shows elevated values before the precipitation event and missing retrievals during rainfall. This behavior is expected because precipitation and cloud cover strongly attenuate lidar signals and prevent reliable direct-sun AOD retrievals by the sun photometer used for cross-validation.

4. b. CIMEL Sunphotometer

Aerosol optical depth (AOD) measured by the CIMEL sun photometer at 1020 nm, the spectral channel closest to the 1064 nm operating wavelength of the ceilometer, is presented in Figure 7 for the two study periods. As noted previously, the CIMEL sun photometer and ceilometer are collocated on the roof of the Grove School of Engineering at The City College of New York (CCNY). The ceilometer-derived and CIMEL-derived AOD values exhibit similar temporal variability; however, the ceilometer-derived AOD is generally lower. This difference is expected because the ceilometer AOD is obtained by integrating aerosol extinction within the planetary boundary layer (PBL), whereas the sun photometer provides a column-integrated AOD measurement from the instrument location to the top of the atmosphere.

Differences between the two retrievals may also arise from their distinct measurement geometries. The sun photometer measures direct solar irradiance along a slant path toward the Sun, while the ceilometer samples the atmosphere along a near-vertical path. Consequently, the instruments may observe different, although spatially proximate, air masses, particularly at greater distances from the measurement site. Nevertheless, their overall temporal patterns are consistent. The lines connecting data points in Figure 7 are included solely as visual guides and do not represent continuous observations. Both instruments indicate higher AOD during June–July 2023 than during June–July 2020. In 2023, AOD frequently exceeded 1.0, whereas values in 2020 only rarely exceeded 0.4. Additional comparisons are provided in the Discussion section.

To investigate the origin and sources of the dense aerosols (such as June 7^th and others), we performed backward trajectory analysis using the HYSPLIT model. The resulting trajectories are shown below. The backward trajectories corresponding to the two spikes observed in the LIDAR-retrieved AOD on June 7 and June 22 (Figure 6) are plotted in Figure 8a and Figure 8b, respectively. From Figure 8a, it is evident that the aerosols within the PBL predominantly originated from Canada, coinciding with particularly severe wildfire conditions during this period (Natural Resources Canada). Figure 8b illustrates clouds formed from aerosol plumes entering the PBL above New York City, within the altitude range of 0.5 km to 1.5 km, originating from outside the local PBL, specifically from northeastern regions of Canada. These aerosols in combination with the marine aerosols present at approximately 1.5 km altitude above the ocean, formed clouds which were transported over the city and precipitated towards the end of the day.

The trajectories clearly indicate that air parcels transported from Canadian regions underwent hygroscopic interactions and heterogeneous mixing with marine aerosols, resulting in cloud formation before being transported back inland as confirmed by the LIDAR vertical profiles illustrated in Figure 8b.

4. c. PM2.5 Measured by TEOM

As mentioned in section 2 of this paper, the PM2.5 data is available in near real-time on the Air Monitoring Website of the New York State Department of Environmental Conservation (NYSDEC). The averaged hour concentrations for the months of June and July of 2020 and 2023 are illustrated in Figure 9. The PM2.5 concentration units are in µg/m³ LC (Local Conditions depending on actual temperature and pressure at the sampling time). Note that the TEOM is located on top of the CCNY Administration building at an altitude of 10 m AGL.

It is already evident in the above figures that the anthropogenic aerosols concentrations are consistently higher in 2023 as compared to 2020. However, more analysis are conducted in the Discussion section.