Atmospheric vertical structure variations during severe aerosol pollution events based on lidar observations

During severe haze events, the boundary layer exhibits a complex vertical structure, while high aerosol loadings hinder high-resolution temperature and humidity measurements. To address this, a Raman–Mie lidar and retrieval algorithms for temperature, humidity, and aerosol optical properties were developed at Xi'an University of Technology, enabling high-resolution profiling of haze vertical structures. A 12 d haze episode was continuously monitored from formation to dissipation, providing detailed spatiotemporal variations of temperature, relative humidity, and aerosols. The boundaries of temperature inversion (TI) and aerosol layers were identified using a threshold method. The results revealed a strong coupling between aerosols and temperature during pollution evolution. Dome and stove effects were observed, with possible coexistence and interaction. The top of a decreasing-type aerosol layer formed a stratified dome-like structure that constrained vertical diffusion, with the temperature gradient of the elevated TI varying inversely with its depth. Both TI strength and humidity were strongly correlated with surface PM_2.5 concentrations. Surface-based TI exhibited a clear diurnal variation, with TI peaks preceding aerosol peaks. The results suggest that strong elevated TI and weak turbulence in the lower layer may facilitate aerosol accumulation. Cloud layers not only suppress radiative heating but may also enhance near-surface humidity through virga processes, which may be conducive to increases in PM_2.5 concentrations. During the dissipation stage, the rapid breakdown of TI and enhanced solar heating were critical for pollutant reduction, while efficient horizontal transport facilitated the complete clearance of aerosols within the boundary layer.

1 Introduction

Air pollution control has long been a key environmental issue accompanying rapid economic development and remains closely linked to public health (Fu et al., 2014; Wang et al., 2023; Song et al., 2017; Di et al., 2021). Aerosols, consisting of fine solid particles and liquid droplets suspended in the atmosphere (Di et al., 2021), represent a key medium for exploring pollution formation and evolution processes. Although excessive aerosol emissions resulting from industrial activities and intensive human socioeconomic processes are the fundamental drivers of frequent pollution events (Zhang et al., 2012, 2013; Le et al., 2020; Wu et al., 2025), meteorological conditions often act as the dominant triggers and regulators of these events (Su et al., 2004; Zhang et al., 2016; Zhu et al., 2018; Zhong et al., 2019; Huang et al., 2021; Jiang et al., 2021; Prasad et al., 2022).

There exists a strong two-way feedback between the vertical distribution of aerosols and meteorological factors (Zhong et al., 2018, 2019). Meteorological conditions unfavorable for aerosol dispersion, such as TI, weak winds, and high relative humidity, are generally considered typical features that exacerbate pollution severity (Zhong et al., 2018, 2019; Miao et al., 2019; Liu et al., 2022; Shao et al., 2023; Zhou et al., 2024). Temperature inversion (TI), characterized by an increase in temperature with height, directly alters atmospheric buoyancy and suppresses vertical mixing (Su et al., 2020). Weak winds restrict horizontal dispersion of aerosols, further facilitating their accumulation (Sekuła et al., 2021). Meanwhile, the hygroscopic growth of aerosols intensifies with increasing relative humidity, serving as a major driver of rapid pollution deterioration (Zhong et al., 2018, 2019). These factors interact synergistically, jointly contributing to the formation and development of surface haze pollution. Huang et al. (2021) investigated the vertical thermal structures under polluted conditions and reported pronounced heating above the planetary boundary layer (PBL), accompanied by cooling near the surface. This pattern, known as the aerosol dome effect (Ding et al., 2016), enhances surface pollutant accumulation by stabilizing the lower atmosphere. Conversely, when aerosols are concentrated near the surface, their absorption of solar radiation heats the lower layer, increases turbulent mixing, and partially alleviates surface pollution – a process termed the aerosol stove effect (Wang et al., 2018). Although numerous studies have explored the interactions between aerosols and meteorological conditions during pollution episodes, most rely on discontinuous or large-scale inferential analyses rather than direct high-resolution observations. Continuous, high spatiotemporal-resolution measurements of boundary-layer structures – particularly temperature and humidity profiles – remain scarce due to observational limitations. This lack of detailed vertical observations hampers a comprehensive understanding of the fine-scale evolution of boundary-layer thermodynamic structures and their coupling with aerosol processes during haze development.

Lidar, characterized by high accuracy and high spatiotemporal resolution, serves as an effective tool for obtaining vertical profiles of temperature, humidity, and aerosols, thereby enabling refined investigations of small-scale atmospheric processes and extreme weather phenomena (Wulfmeyer et al., 2015; Girolamo et al., 2016; Lange et al., 2019; Yi et al., 2021; Li et al., 2022, 2023, 2025). Under clear-sky conditions, lidar can provide highly precise temperature and humidity measurements. However, under high aerosol loading (i.e., during haze pollution), signal crosstalk may induce substantial deviations in retrieved temperature profiles (Su et al., 2013; Li et al., 2022, 2023, 2025). This issue mainly arises because strong Mie scattering signals contaminate the rotational Raman channel used for temperature detection, leading to retrieval errors. In addition, discrepancies in the geometric overlap factors among Raman channels can introduce systematic biases in temperature and humidity retrievals. These limitations make accurate thermodynamic profiling particularly challenging in polluted atmospheres. To address these challenges, Xi'an University of Technology (XUT) developed a multi-wavelength Raman–Mie scattering lidar system. Di, Li, and co-workers conducted comprehensive studies on temperature and humidity retrieval techniques under hazy and cloudy conditions and proposed correction algorithms for rotational Raman temperature retrievals and geometric overlap effects. These advancements enabled high-precision retrievals of temperature, humidity, and aerosol profiles within clouds and throughout the atmospheric boundary layer (Li et al., 2022, 2023, 2025).

In the winter of 2023, a severe haze pollution episode was continuously observed using the XUT Raman–Mie scattering lidar system. By applying the aforementioned retrieval algorithms, high-accuracy profiles of temperature, humidity, and aerosol optical properties within the haze layer were obtained, allowing detailed characterization of the vertical evolution of the atmosphere during the formation, development, and dissipation stages of the event. Multiple auxiliary datasets, including radiosonde soundings and ERA5 reanalysis data, were integrated to examine the coupling between meteorological conditions and pollution evolution. The observations revealed a distinct dome-shaped stratification of aerosols induced by solar radiation, demonstrating the coexistence of both the dome effect and the stove effect. Furthermore, the temporal evolution of TIs and their regulatory influence on the vertical distribution of aerosols were clearly captured throughout the pollution episode. Correlation analyses were conducted between surface PM_2.5 concentrations and thermodynamic variables, such as temperature and relative humidity, alongside assessments of the diurnal variations in temperature gradients.

2 Data, observation locations, and methods

2.1 Data

The primary datasets employed in this study include lidar-derived profiles of temperature, humidity, and aerosol optical properties; radiosonde soundings; surface meteorological observations; surface particulate matter concentrations (PM₁₀ and PM_2.5) obtained from air quality monitoring stations; surface radiation measurements; and vertical atmospheric profiles from the ERA5 reanalysis dataset.

The Raman–Mie scattering lidar employed in this study was developed by Xi'an University of Technology in 2021. The system uses a 354.7 nm Nd:YAG laser as the excitation source and a 400 mm-aperture Cassegrain telescope to collect atmospheric backscatter signals. The complete overlap blind zone of the system is approximately 120 m. It provides a minimum vertical resolution of 3.75 m and a temporal resolution of 2 min. Atmospheric temperature profiles are retrieved using a rotational Raman method, whereas water vapor profiles are derived via a combined vibrational-rotational Raman approach. A vibrational Raman channel for nitrogen further enables accurate retrievals of aerosol extinction and backscatter coefficients using the Raman technique, with retrieval uncertainties below 15 %. For temperature retrievals, a 30 min temporal accumulation and 60-point smoothing are applied to enhance the signal-to-noise ratio, while water vapor retrievals use a 10 min accumulation and 30-point smoothing. Relative humidity is calculated from the lidar-derived temperature and water vapor profiles. Comparisons with co-located radiosonde measurements show that the temperature and water vapor uncertainties are within 1 K and 0.5 g kg⁻¹, respectively. The main technical specifications of the lidar system are summarised in Table 1.

Table 1 Main technical specifications.

HFOV, half field of view; TA, telescope aperture; EW, excitation wavelength; EPP, energy per pulse; PD, pulse duration; FWHM, full width at half maximum; ESSR, elastic scattering suppression ratio.

Radiosonde data were obtained from the Xi'an Jinghe Meteorological Basic Station, where soundings were launched daily at 07:15 and 19:15 Beijing Time. The radiosonde model used was GTS13, transmitting measurements at a frequency of 1 Hz, with a vertical resolution of 6 m in the lower atmosphere and a maximum detection height of 30 km. When the ambient temperature exceeds −80 °C, the temperature measurement accuracy is ±0.2 °C. The relative humidity measurement error is less than 6 % under ambient temperatures above −20 °C, and for atmospheric pressures above 500 hPa, the pressure measurement error is below 1.5 hPa. The instrument also provides horizontal wind speed and direction.

Surface observations were collected from ground-based meteorological stations and particulate matter (PM) monitors, providing hourly measurements of near-surface temperature, humidity, atmospheric pressure, visibility, and PM concentrations.

PM₁₀ and PM_2.5 concentrations were measured using β-ray absorption particulate matter monitors (LGH-01B and LGH-01E). Both instruments operate over a concentration range of 0–1000 µg m⁻³, with a minimum sampling interval of 30 min and a sampling flow rate of 16.7 L min⁻¹. The lower detection limit is $<\mathrm{2}\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$, and the flow stability is within 2 %. Surface meteorological variables were obtained from an automatic weather station (DZZ5). Atmospheric pressure is measured over a range of 500–1100 hPa with a resolution of 0.1 hPa and a maximum permissible error of ±0.3 hPa. Air temperature is measured over a range of −50 to 50 °C, with a resolution of 0.1 °C and a maximum permissible error of ±0.2 °C. Relative humidity is measured over 5 %–100 %, with a resolution of 1 %, and maximum permissible errors of ±3 % for RH≦80 % and ±5 % for RH>80 %.Surface solar irradiance was measured using a TBQ-2-UMB pyranometer, which has a spectral response range of 0.3–3 µm and a measurement range of 0–2000 W m⁻². The hourly measurement uncertainty is approximately 8 %. These surface measurements were co-located with the radiosonde launch site to ensure data consistency.

ERA5 is the fifth-generation global atmospheric reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). It provides hourly estimates of a wide range of atmospheric, land, and oceanic variables on a regular latitude–longitude grid with a horizontal resolution of approximately 31 km (0.25°×0.25°) and 137 model levels from the surface up to 0.01 hPa. ERA5 employs a four-dimensional variational data assimilation system to combine model output with a large number of quality-controlled observations, generating a physically consistent, spatially complete, and temporally continuous representation of the Earth system. The dataset includes variables such as temperature, humidity, horizontal and vertical winds, surface pressure, precipitation, and radiation, available at the surface and multiple pressure or model levels. ERA5 data were obtained from the Copernicus Climate Data Store (CDS, https://cds.climate.copernicus.eu, last access: 1 July 2025). In this study, we focused on key meteorological parameters, including isobaric maps, horizontal wind vectors, and vertical motion, which were used to characterize the background atmospheric conditions.

2.2 Observation locations

The above instruments were located at the Jinghe National Meteorological Station in Xi'an, Shaanxi Province, China (34°26^′ N, 108°58^′ E), situated on the north bank of the Wei River in the north-central region of the Guanzhong Plain (33°42^′–34°45^′ N, 107°40^′–109°49^′ E, about 400 m above sea level). The geographical location of the study area is shown in Fig. 1. The Guanzhong Plain is a typical elongated basin-oriented east–west, bordered by the Qinling Mountains to the south and the northern Loess Plateau to the north. Elevations are higher along the southern and northern margins and gradually decrease toward the east. This enclosed or semi-enclosed topography restricts both horizontal and vertical air exchange, impeding the dispersion of atmospheric pollutants. In winter, the influence of the westerlies weakens, and the frequency and intensity of cold air outbreaks decrease. Under these conditions, stable stratification and TIs occur more frequently, and the boundary layer is typically lower, further limiting pollutant dilution and transport. Consequently, the Guanzhong Plain is prone to prolonged, large-scale, and severe air pollution episodes during winter.

Figure 1 Geographical coverage around observation site (31°30^′–40°N, 105°–112°E). The red dot indicates the location of the Jinghe National Meteorological Station in Xi'an (34°26^′ N, 108°58^′ E)

2.3 Lidar data retrieval algorithm

2.3.1 Retrieval algorithms for temperature, humidity, and backscatter coefficient within haze layers

The fundamental principle for probing atmospheric vertical structures using a Raman–Mie scattering lidar is based on the lidar equation, which involves the theories of Mie–Rayleigh scattering, rotational Raman scattering, and vibrational Raman scattering. This can be represented by Eqs. (1)–(3). Specifically, rotational Raman scattering is used for temperature measurements, while Mie–Rayleigh scattering and vibrational Raman scattering are employed for retrieving aerosol optical properties and water vapor profiles.

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\begin{array}{rl}{\displaystyle }{P}_{\mathrm{r}}(T,z,{\mathit{\lambda }}_{\mathrm{r}})=& {\displaystyle \frac{C_{\mathrm{r}}}{z^{\mathrm{2}}}}{\mathit{\beta }}_{\mathrm{r}}(T,z,{\mathit{\lambda }}_{\mathrm{r}})\\ {\displaystyle }& {\displaystyle }\times \exp \left\{-{\int }_{\mathrm{0}}^z\left[\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\mathrm{e}})+\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\mathrm{r}})\right]\mathrm{d}{z}^{\prime }\right\}\end{array}\end{array}$$

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\begin{array}{rl}{\displaystyle }{P}_{\mathrm{v}}(z,{\mathit{\lambda }}_{\mathrm{v}})=& {\displaystyle \frac{C_{\mathrm{v}}}{z^{\mathrm{2}}}}{\mathit{\beta }}_{\mathrm{v}}\left(z,{\mathit{\lambda }}_{\mathrm{v}}\right)\\ {\displaystyle }& {\displaystyle }\times \exp \left\{-{\int }_{\mathrm{0}}^z\left[\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\mathrm{e}})+\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\mathrm{v}})\right]\mathrm{d}{z}^{\prime }\right\}\end{array}\end{array}$$

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\begin{array}{rl}{\displaystyle }{P}_{\mathrm{e}}(z,{\mathit{\lambda }}_{\mathrm{e}})=& {\displaystyle \frac{C_{\mathrm{e}}}{z^{\mathrm{2}}}}\left({\mathit{\beta }}_{\mathrm{a}}\right(z,{\mathit{\lambda }}_{\mathrm{e}})+{\mathit{\beta }}_{\mathrm{m}}(z,{\mathit{\lambda }}_{\mathrm{e}}\left)\right)\\ {\displaystyle }& {\displaystyle }\times \exp \left\{-\mathrm{2}{\int }_{\mathrm{0}}^z\left[\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\mathrm{e}})\right]\mathrm{d}{z}^{\prime }\right\}\end{array}\end{array}$$

where P_r, P_v, and P_e denote the rotational Raman, vibrational Raman, and elastic backscattering signals, respectively. C represents the system constant of the lidar, which includes the laser power, telescope receiving area, optical transmittance, and detector quantum efficiency. λ, T, and z refer to the wavelength, temperature, and height, respectively, while α and β denote the volume extinction and backscattering coefficients. The subscripts e, rh, rl, vw, and vn correspond to the Mie–Rayleigh signal, high-quantum-number rotational Raman signal, low-quantum-number rotational Raman signal, water vapor vibrational Raman signal, and nitrogen vibrational Raman signal, respectively, whereas m and a represent the molecular and aerosol components. The atmospheric extinction coefficient can be derived from the nitrogen vibrational Raman signal (P_vn), while the backscatter ratio (BR) is obtained in combination with the Mie–Rayleigh scattering signal. Finally, the aerosol backscatter coefficient can be retrieved, as expressed in Eqs. (4)–(6).

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{\displaystyle }{\displaystyle }{\mathit{\alpha }}_{\mathrm{a}}(z,{\mathit{\lambda }}_{\mathrm{e}})={\displaystyle \frac{\mathrm{d}\left[\ln \left(\frac{N\left(z\right)}{z^{\mathrm{2}}{P}_{\text{vn}}}\right)\right]/\mathrm{d}z-{\mathit{\alpha }}_{\mathrm{m}}(z,{\mathit{\lambda }}_{\mathrm{e}})-{\mathit{\alpha }}_{\mathrm{m}}(z,{\mathit{\lambda }}_{\text{vn}})}{\mathrm{1}+({\mathit{\lambda }}_{\mathrm{e}}/{\mathit{\lambda }}_{\text{vn}}{)}^K}}\text{(5)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }\text{Br}(z,{\mathit{\lambda }}_{\mathrm{e}})=& {\displaystyle }F{\displaystyle \frac{C_{\text{vn}}}{C_{\mathrm{e}}}}{\displaystyle \frac{P_{\mathrm{e}}(z,{\mathit{\lambda }}_{\mathrm{e}})}{P_{\text{vn}}(z,{\mathit{\lambda }}_{\text{vn}})}}\\ {\displaystyle }& {\displaystyle }\times \exp \left\{{\int }_{\mathrm{0}}^z\left[\mathit{\alpha }\right(z^{\prime },{\mathit{\lambda }}_{\mathrm{e}})-\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\text{vn}}\left)\right]\mathrm{d}{z}^{\prime }\right\}\end{array}\text{(6)}& {\displaystyle }& {\displaystyle }{\displaystyle }{\displaystyle }{\mathit{\beta }}_{\mathrm{a}}(z,{\mathit{\lambda }}_{\mathrm{e}})=\text{Br}(z,{\mathit{\lambda }}_{\mathrm{e}})\cdot {\mathit{\beta }}_{\mathrm{m}}(z,{\mathit{\lambda }}_{\mathrm{e}})-{\mathit{\beta }}_{\mathrm{m}}(z,{\mathit{\lambda }}_{\mathrm{e}})\end{array}$$

In Eq. (4), N(z) and K represent the atmospheric molecular number density and the Ångström exponent, respectively, which is typically taken as 1. In Eq. (5), F denotes the system constant used for calculating BR. The rotational Raman temperature retrieval is based on the correlation between the rotational Raman signal and atmospheric temperature. By extracting two rotational Raman backscattering signals with different temperature sensitivities, a functional relationship between their signal ratio and temperature can be established to derive the atmospheric temperature profile. However, the rotational Raman channels are highly susceptible to elastic scattering cross-talk, which may introduce biases in the temperature retrieval. In Eq. (7), A and B are system constants used for temperature retrieval. M_h and M_l represent the elastic scattering cross-talk coefficients for the high- and low-quantum-number rotational Raman channels, respectively. The water vapor and nitrogen vibrational Raman scattering signals are extracted, and their ratio is calculated. Based on radiosonde data obtained under the same temporal and spatial conditions, a regression relationship for data retrieval is constructed, enabling the lidar measurement of the water vapor mixing ratio and, consequently, the derivation of relative humidity, as expressed in Eq. (8). The aforementioned measurement principles are well established, and detailed descriptions can be found in previous studies (Ansmann et al., 1992; Li et al., 2023, 2025; Chazette et al., 2025).

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }T\left(z\right)=\\ {\displaystyle }& {\displaystyle }\mathrm{A}/\left\{\ln \left[{\displaystyle \frac{C_{\text{rl}}\left(z\right)}{C_{\text{rh}}\left(z\right)}}{\displaystyle \frac{P_{\text{rh}}(T,z,{\mathit{\lambda }}_{\text{rh}})+M_{\mathrm{h}}\cdot P_{\mathrm{e}}(z,{\mathit{\lambda }}_{\mathrm{e}})}{P_{\text{rl}}(T,z,{\mathit{\lambda }}_{\text{rl}})+M_{\mathrm{l}}\cdot P_{\mathrm{e}}(z,{\mathit{\lambda }}_{\mathrm{e}})}}\right]-\mathrm{B}\right\}\end{array}\text{(8)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }W\left(z\right)=& {\displaystyle \frac{C_{\text{vn}}}{C_{\text{vw}}}}{\displaystyle \frac{P_{\text{vw}}(z,{\mathit{\lambda }}_{\text{vw}})}{P_{\text{vn}}(z,{\mathit{\lambda }}_{\text{vn}})}}\\ {\displaystyle }& {\displaystyle }\times \exp \left\{{\int }_{\mathrm{0}}^z\left[\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\text{vw}})-\mathit{\alpha }(z^{\prime },{\mathit{\lambda }}_{\text{vn}})\right]\mathrm{d}{z}^{\prime }\right\}\end{array}\end{array}$$

Under atmosphere conditions with cloud layers or high aerosol concentrations, due to the close spectral intervals, strong elastic scattering can cause elastic scattering cross-talk in the low-quantum-number rotational Raman channel, and even in the high-quantum-number channel. This means that a portion of the elastic scattering signal is mixed into the rotational Raman channels, which leads to significant deviations in data retrieval, and in severe cases, preventing accurate temperature retrieval within the haze layer. Since haze layers are generally located at low height, the geometric overlap factors of each lidar channel are not identical, which also causes large uncertainties in temperature retrievals in the lower atmospheric region. Therefore, temperature measurements in haze layers mainly include three processes: geometric overlap factor correction, rotational Raman ratio correction, and temperature retrieval.

The atmospheric temperature correction technique is based on establishing a linear relationship between the backscatter ratio and the elastic-scattering crosstalk ratio, which allows the rotational Raman ratio to be corrected using the measured backscatter ratio. Figure 2 illustrates the temperature retrieval process based on this algorithm. First, a high-altitude region unaffected by the geometric overlap factor is selected, and radiosonde temperature profiles are used for system calibration. Under clear-sky and dry near-surface conditions, the theoretical rotational Raman ratio is derived from radiosonde data (black dash–dotted line in Fig. 2a). By comparing it with the measured Raman ratio (blue solid line in Fig. 2a), the geometric overlap factor is obtained (red solid line in Fig. 2b), enabling overlap correction. Next, the theoretical Raman ratio under strong elastic-scattering conditions is derived from radiosonde data (black dash–dotted line in Fig. 2d). The measured Raman ratio (black solid line in Fig. 2d) is then used to calculate the elastic-scattering crosstalk ratio. A linear regression between the backscatter ratio and the crosstalk ratio is performed to determine the system calibration constant (Fig. 2e). Using this constant together with the measured backscatter ratio, the rotational Raman ratio is corrected, enabling retrieval of the atmospheric temperature profile within the elastic-scattering region (Fig. 2f). A linear regression between the corrected lidar temperature and the radiosonde temperature yields a coefficient of determination (R²) of approximately 0.8 (Fig. 2g), indicating good agreement between the two measurements. A detailed description of the algorithm can be found in Li et al. (2025).

Figure 2 Atmospheric temperature correction. (a) Measured (blue solid) and theoretical (black dash–dotted) rotational Raman ratios. (b) Overlap-related quantities, including the measured ratio (black dashed) and the overlap function (red solid). (c) Range-square-corrected signals from the elastic (Mie–Rayleigh; thick red solid), nitrogen vibrational Raman (thin purple solid), and high- and low-quantum-number rotational Raman (black dash–dotted and blue dashed) channels. (d) Backscatter ratio (blue dash–dotted) and rotational Raman ratios (measured ratio, black solid; theoretical ratio, black dash–dotted). (e) Linear regression between the backscatter ratio (BR) and the cross-talk ratio. (f) Backscatter coefficient (BC) and temperature profiles derived from lidar measurements and radiosonde observations. The black dash–dotted line denotes the radiosonde temperature, while the red solid and blue dashed lines represent the corrected and uncorrected lidar temperature profiles, respectively. Shaded areas indicate the corresponding uncertainties. (g) Linear regression between radiosonde and lidar temperatures.

2.3.2 Retrieval algorithms for vertical velocity and buoyancy acceleration

The vertical velocity used in this study was derived from the vertical pressure tendency provided by ERA5. Atmospheric static stability is governed by buoyancy, which in turn influences the vertical distribution of aerosols. In this study, the buoyancy acceleration was estimated using an approach based on the environmental virtual potential temperature. This method assumes small perturbations and linearization, making it suitable for weakly disturbed conditions. Under strong perturbations, however, it may underestimate the actual buoyancy. Nevertheless, the estimated values can reasonably capture the overall trend of changes in thermodynamic stability. Buoyancy is expressed as (Su et al., 2020)

$$\begin{array}{}\text{(9)}& B={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{z}_{\mathrm{0}}}{\mathrm{d}{t}^{\mathrm{2}}}}=-g\mathrm{\Delta }z{\displaystyle \frac{\mathrm{1}}{\mathit{\theta }}}{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}{z}_{\mathrm{0}}}}\end{array}$$

where z₀ denotes the height of the air parcel, and t represents time. θ is the virtual potential temperature of the environment, which is calculated using radiosonde pressure combined with the lidar-measured atmospheric temperature and water vapor mixing ratio. A positive buoyancy indicates a convective atmospheric layer, whereas a negative buoyancy corresponds to a stable layer. When the buoyancy approaches zero, the atmospheric column is neutral.

According to the theory of error propagation, the uncertainty of buoyancy acceleration can be expressed as:

$$\begin{array}{}\text{(10)}& {\mathit{\delta }}_B=\sqrt{{\left({\displaystyle \frac{\partial B}{\partial \mathrm{\Gamma }}}\right)}^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{\Gamma }}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial B}{\partial \mathrm{\Gamma }}}\right)}^{\mathrm{2}}{\mathit{\delta }}_{\mathit{\theta }}^{\mathrm{2}}}\end{array}$$

In the formula, $\mathrm{\Gamma }=\mathrm{d}\mathit{\theta }/\mathrm{d}{z}_{\mathrm{0}}$. The lidar system exhibits a high signal-to-noise ratio within the detection range below 4 km at night (2 km during daytime), and the temperature uncertainty is typically less than 0.5 K. Selecting typical boundary layer conditions with a virtual temperature of approximately 280 K and a conservative uncertainty of 0.5 K, when the vertical resolution is 37.5 m, the uncertainty of the temperature gradient obtained is approximately 0.019 K m⁻¹, and the uncertainty of the buoyancy acceleration is calculated to be approximately 0.025 m s⁻². This value represents the worst-case scenario assuming completely uncorrelated errors. Lidar temperature errors exhibit strong vertical correlation, and buoyancy is essentially derived from vertical gradients. Therefore, the effective uncertainty in buoyancy is substantially smaller. Although this uncertainty may affect very weak buoyancy signals near neutral conditions, the identification of stable layers and strong inversions is robust.

2.3.3 Method for Determining Aerosol Boundary and TI Characteristic

The aerosol layer boundary was identified using a gradient-based method applied to BR. The fundamental principle of this method is to locate positions exhibiting pronounced discontinuities in the BR profile by calculating its vertical derivative. As this study primarily focuses on the relationship between the aerosol vertical structure and temperature distribution, only the upper boundary of the aerosol layer (hereafter referred to as UABL) was determined, without further distinction of the stable boundary layer, mixed layer, and residual layer.

First, the first derivative of the BR profile (i.e., BR gradient) was calculated, and the negative-gradient intervals of BR (hereafter referred to as NIBRs) were identified based on the zero-crossing points of the derivative. To eliminate local perturbations and high-frequency noise, threshold conditions (>60 m) were applied to both the NIBR lengths and the spacing between adjacent NIBRs, ensuring that only valid NIBRs were retained. Similarly, the positive-gradient intervals of BR (hereafter referred to as PIBRs) were identified through zero-crossing detection and filtered with the same threshold criteria. Second, gradient thresholds were defined for both PIBRs ($>\mathrm{3}{\mathrm{km}}^{-\mathrm{1}}$) and NIBRs ($<-\mathrm{3}{\mathrm{km}}^{-\mathrm{1}}$). In addition, a minimum BR deviation threshold within each interval ($\left|\mathrm{\Delta }\text{Br}\right|>\mathrm{0.35}$, defined as the absolute deviation between the maximum and minimum BR within the interval) was imposed to exclude pseudo-transition zones with weak gradients. When the length of a positive interval satisfying the gradient threshold exceeded 30 m, it indicated a notable enhancement of aerosol concentrations within that layer. In such cases, the adjacent negative interval immediately above was retained even if it did not meet the gradient threshold. Finally, the UABL was determined within NIBRs using a BR threshold of 1.3. Starting from the top of each NIBR, one-eighth of its length was used as the reference range. If the maximum BR within this range exceeded 1.3, the height corresponding to one-eighth of the NIBR length below the top was defined as the UABL. Otherwise, the height at which the BR first exceeded 1.3 within NIBRs was identified as the UABL.

Figure 3 Extraction of UABL and TI characteristics. Panels (a)–(d) show (a) BR, (b) BR gradient, (c) temperature, and (d) temperature gradient. In panels (a)–(d), the black curves represent the BR, BR gradient, temperature, and temperature gradient, respectively. Asterisks indicate the identified aerosol layer boundaries. In panel (b), the blue solid and red dashed shaded regions denote NIBR and PIBR, respectively. In panels (c) and (d), the red dashed shaded regions highlight the TI layers.

TI features were identified using a gradient-based method. Unlike conventional gradient methods, the raw temperature profile was first subjected to linear feature fitting prior to calculating the vertical temperature gradient. The fitting employed the piecewise linear interpolation method proposed by Fochesatto (2015), which constructs an interpolation function with variable segment lengths and minimizes an associated error function to obtain optimal fits for each segment, thereby accurately resolving TI layers.

The first derivative of the fitted temperature profile (i.e., the vertical temperature gradient) was then calculated, and the positive-gradient intervals of temperature (hereafter referred to as PITs) along with their adjacent gap layers were identified using zero-crossing points. Threshold conditions (>60 m) were applied to both PIT lengths and gap layers to retain only valid segments. Gap layers shorter than the threshold were merged with adjacent PITs to form a single TI layer. Conversely, PITs shorter than the threshold, when adjacent to gap layers exceeding the threshold, were discarded as likely resulting from minor perturbations or weak TI features. To prevent underestimation of TI thickness, additional criteria were applied to isolated PITs. If a gap layer exceeded the threshold length and its mean vertical gradient exceeded $-\mathrm{0.5}\mathrm{K}{\mathrm{km}}^{-\mathrm{1}}$, or the maximum temperature deviation within the gap layer was less than 0.5 K, indicating nearly uniform temperature with height, the adjacent PITs were merged. PITs satisfying these conditions were thus considered valid TI layers. The base and top of each TI layer were defined as the TI base and TI top, respectively. The vertical distance between TI base and top represents TI thickness. The equivalent vertical temperature gradient was then calculated from the temperature difference within TI layer.

3 Results and discussion

3.1 Overview of the haze pollution event

A prolonged haze pollution episode occurred in Xi'an in late December 2023, lasting for 12 d. Figure 4 illustrates the surface PM concentrations and meteorological parameters during this persistent pollution event. As shown in the figure, the surface aerosol concentrations began to accumulate on 22 December and continued to increase, reaching a peak in PM_2.5 concentrations on the evening of 28 December, which persisted until noon on 30 December. A brief reduction of pollutants occurred in the afternoon of 30 December, followed by a rapid re-accumulation during the night. Around noon on 31 December, the surface PM_2.5 concentrations decreased to approximately 100 µg m⁻³ and remained at this level until noon on 2 January 2024, after which the pollution entered its dissipation stage, marking the end of the event. According to the temporal variation of PM concentrations, this pollution process can be divided into three aerosol accumulation stages (AASs), as indicated by the red dashed boxes in Fig. 4a, with the first two stages representing the dominant accumulation phases. As shown in Fig. 4b, surface temperature and relative humidity exhibited regular diurnal fluctuations beginning on 22 December and showed an overall increasing trend. However, from noon on 28 December to noon on 29 December, both temperature and humidity remained nearly constant. Starting on 30 December, they resumed their diurnal variations, while humidity showed a gradual decreasing trend. The surface pressure showed an inverse relationship with the aerosol concentrations, reaching its minimum when pollution was most severe. Figure 4c shows that horizontal visibility gradually decreased from the evening of 22 December to noon on 30 December, with the lowest visibility below 1 km. From the evening of 31 December to noon on 2 January 2024, visibility again showed a decreasing trend, reaching a minimum of about 5 km. As the pollution dissipated, horizontal visibility improved rapidly.

Figure 4 Surface PM and meteorological parameters. Panels (a–c) show (a) PM concentrations, (b) surface meteorological parameters, and (c) visibility. In panel (a), the black and blue dotted curves represent the PM_2.5 and PM₁₀ concentrations, respectively. In panel (b), the black, blue dotted, and pink dashed curves represent the temperature, relative humidity, and pressure, respectively.

Figure 5 illustrates the synoptic circulation pattern at 900 hPa over the study region (101°–115° E, 28°–40° N) during the haze pollution episode. During this period, the Guanzhong Plain was dominated by a weak high-pressure system, characterized by a gentle pressure gradient, stagnant moisture transport, and weak wind variability. These conditions resulted in poor horizontal and vertical ventilation, thereby inhibiting the dispersion of pollutants. The persistent weak synoptic forcing and stable stratification provided favorable meteorological conditions for aerosol accumulation over the region. To assess the role of large-scale transport in the development of the pollution episode, 48 h backward air-mass trajectories at multiple altitudes (0.5, 1.5, and 3.0 km) were analyzed for 28 December, which was characterized by severe pollution. The results indicate that airflows at approximately 0.5 and 1.5 km were predominantly advected from southern regions toward the observation site, with relatively short transport pathways and limited residence time over major upwind pollution source areas. In contrast, air masses at 3.0 km mainly originated from the west to northwest and exhibited weak coupling with the near-surface aerosol layer. Consequently, large-scale regional transport is unlikely to be the dominant factor governing the near-surface aerosol loading during this period. Moreover, strict emission control measures were implemented in Xi'an, including traffic restrictions and constraints on industrial emissions, which may have further reduced local primary emissions.

Figure 5 900 hPa geopotential height (black contours, gpm), temperature (red contours, °C), and wind vectors at 00:00 UTC on (a) 24 December, (b) 26 December, (c) 28 December, and (d) 30 December 2023. The blue triangle marks the location of Xi'an.

3.2 High-resolution vertical structure of the pollution episode

The evolution of this pollution event was observed by both ground-based meteorological instruments and a lidar system. The lidar provided vertical profiles of aerosols and clouds from 0 to 10 km, as well as temperature and relative humidity from 0 to 4 km, with profiles recorded every 2 min (in Fig. 6). Periods of lidar maintenance occurred during 23–24, 26, 27 December 2023, and 2 January 2025, during which the system was not operating, resulting in no data. In Fig. 6c–e, the white areas in the upper layers represent regions of invalid lidar detection, primarily due to severe low-level haze reducing the lidar penetration capability and thus lowering the signal-to-noise ratio of the Raman signal. Figure 6a and c show the 1064 nm Mie–Rayleigh backscatter signal and the 355 nm backscatter coefficient, respectively, with Fig. 6c also indicating the top height of the aerosol layer. Figure 6b presents net and total radiation fluxes observed by ground-based instruments. Figure 6f–h display radiosonde data collected every 12 h.

Figure 6 Spatiotemporal distributions of atmospheric vertical structure observed by lidar and radiosonde: (a) range-square-corrected signal (RSCS) at 1064 nm; (b) radiant exposure; (c) backscatter coefficient at 355 nm; (d) temperature measured by lidar; (e) relative humidity measured by lidar; (f) temperature from radiosonde; (g) relative humidity from radiosonde; (h) wind vectors and wind speed from radiosonde.

As shown in Fig. 6a, three distinct cloud events occurred in the upper atmosphere during the pollution episode, specifically on 22–23, 28–29, and 31 December 2023–1 January 2024. All three events were observed to coincide with reductions in surface radiation. The second and third cloud events were accompanied by virga, during which increased near-surface relative humidity was observed, temporally coinciding with periods of elevated surface pollution. The continuous vertical profile of the lidar Mie–Rayleigh backscatter signal served as the primary basis for the preliminary identification of cloud layers and virga, which could be further corroborated by depolarization ratio measurements or radial velocity observations from millimeter-wave radar (Yi et al., 2021; Zou et al., 2024; Jimenez et al., 2025). Figure 6c clearly shows that the near-surface backscatter coefficient began to increase gradually on the evening of 22 December, reaching high values during 28–30 December, consistent with variations in surface PM concentrations. A temporary decrease occurred on the afternoon of 30 December, followed by a rapid increase in the near-surface backscatter coefficient. From the evening of 31 December to noon on 2 January, another increasing trend was observed, although the values were much lower than before, and by the afternoon of 2 January, the values had decreased to levels comparable to pre-pollution conditions. We divide the pollution episode into two stages: the pollution development stage (21 December to noon on 30 December) and the late pollution dissipation stage (afternoon of 30 December–3 January 2024).

Figure 7 presents the continuous distributions of potential temperature, buoyancy acceleration, and vertical velocity. The vertical velocity used in this study was derived from the vertical pressure tendency provided by ERA5. It is noteworthy that ERA5-derived vertical velocity is not used to quantitatively interpret the fine-scale aerosol stratification observed by the lidar. Instead, it is employed solely to characterize the large-scale synoptic background, such as periods of weak subsidence or ascent that provide a favorable or unfavorable environment for boundary-layer development.

Figure 7 Spatiotemporal distributions of potential temperature, buoyancy acceleration, and vertical velocity: (a) potential temperature; (b) buoyancy acceleration; (c) vertical velocity.

3.2.1 Pollution development stage

The dominant phase of the pollution development stage occurred from 22–28 December. Figure 8 shows the vertical distributions of aerosols and temperature during this period, including profiles in the morning and evening, and daily averages. Over the course of the stage, the vertical distributions of aerosols and temperature underwent persistent changes. At the early stage of pollution, the aerosol vertical structure was primarily characterized by a decreasing profile with height, or alternating patterns of decreasing and well-mixed near-surface layers. The temperature vertical structure exhibited pronounced diurnal variations. During periods of severe pollution, the aerosol vertical structure became well-mixed, and the diurnal variations in both aerosol and temperature distributions were not significant. On the evening of 22 December (in Fig. 8a), the backscatter coefficient below approximately 0.7 km was significantly higher than in the morning, while above 0.7 km it was lower than the morning values. The base heights of the elevated TI layer in the morning, noon, and evening were approximately 1.3, 1.1, and 0.7 km, respectively. Subsequently, a surface-based TI layer with a thickness of about 1 km developed, strongly inhibiting the vertical diffusion of aerosols. It should be noted that the description of the surface-based TI layer here refers only to the region above the reliable temperature retrieval height (>120 m). This surface-based TI persisted until 24 December, resulting in a stable and well-mixed aerosol layer below 0.5 km. On 25 December, surface heating due to solar radiation promoted vertical aerosol mixing, but also led to the formation of a dome-like aerosol layer around 1.1 km, producing a dome effect that enhanced the TI at this height and restricted aerosols below approximately 1.0 km. Notably, under heavy pollution conditions, the lower atmosphere did not experience a strong TI, but the temperature gradient below 0.8 km was about 3.5 K km⁻¹, indicating a stable atmospheric state. The vertical diffusion of aerosols was further limited by the elevated TI layer above, confining aerosols to below the TI height.

Figure 8 Vertical profiles of aerosol backscatter coefficient and temperature from 22–28 December. Panels (a)–(g) show the aerosol backscatter coefficient profiles. The red dashed and blue dash–dotted lines denote the profiles at 07:00 and 19:00 CST, respectively. The black solid line represents the daily mean profile, while the grey solid lines indicate hourly profiles. Red circles and blue squares mark the upper boundary of the aerosol layer in the morning and evening, respectively. Panels (h)–(n) show the temperature profiles with the same plotting conventions. Red circles and blue squares indicate the top heights of the TI in the morning and evening, respectively. Panel (o) shows the daily mean PM_2.5 concentrations.

Figure 9 presents the hourly vertical profiles of aerosols and temperature during the early stage of aerosol accumulation. The UABL height exhibited a strong positive correlation with the TI layer and a negative correlation with surface PM_2.5 concentrations. On 21 December, during the pollution brewing stage, a pronounced stratification was observed at heights of 2–3 km. From the morning of 21 to 22 December, a significant TI layer occurred between approximately 1.5 and 2.5 km, inhibiting the upward transport of aerosols. However, because the TI layer was relatively elevated, vertical transport below its base remained possible, resulting in a decreasing aerosol profile, with higher concentrations near the surface that decreased with height. The average PM_2.5 concentrations near the surface remained below 50 µg m⁻³, indicating good air quality. During the daytime of 22 December, cloud cover significantly reduced downward radiation fluxes (in Fig. 6b), which directly affected near-surface radiative heating and enhanced the stability of the lower atmosphere, limiting the vertical diffusion of water vapor and aerosols. In addition, the residual layer produced a dome effect, causing pronounced warming above approximately 1 km around noon and strengthening the TI in this layer. The elevated TI induced downward buoyancy acceleration, further reducing the UABL height and confining free vertical convection to below approximately 0.9 km (in Fig. 7). Moreover, from 21–22 December, winds below 500 m were predominantly northeasterly at 5–10 m s⁻¹. In winter in Xi'an, such northeasterly winds are typically associated with pollution transport, initiating the first aerosol accumulation on the afternoon of 22 December. After sunset, a stable layer developed below approximately 1 km, favoring the accumulation of aerosols in the near-surface layer. The vertical structure of near-surface aerosols transitioned from a decreasing profile to a well-mixed profile, in which concentrations remained nearly constant with height. On 23 December, weak and variable winds ($<\mathrm{3}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) prevailed near the surface, while cloud cover continued to inhibit radiative transfer, maintaining a stable TI structure. Surface PM_2.5 concentrations continued to increase. Considering the buoyancy conditions during this period, the vertical convection range was substantially compressed. Thus, the subsidence of the TI at the UABL, combined with cloud radiative forcing, appears to jointly constrain the vertical convective scale during the early stage of pollution and may be a primary factor contributing to aerosol accumulation.

Figure 9 The early stage of aerosol accumulation. Panels (a) and (b) show the backscatter coefficient and temperature, respectively. In panel (a), black solid lines represent the backscatter coefficient, and dashed lines indicate the UABL layers, with shaded dashed boxes highlighting the UABL heights. When aerosol stratification occurs, red upward-pointing triangles and black downward-pointing triangles correspond to the upper and lower UABLs, respectively. In panel (b), the blue solid line represents the temperature, while red dashed lines and black solid lines correspond to lidar measurements and radiosonde observations, respectively. Shaded dashed boxes indicate the locations of TI layers. The yellow shaded regions represent the magnitudes of PM_2.5 concentrations.

During the mid-stage of pollution development, the vertical structure of aerosols exhibited a well-mixed profile, although the boundary layer remained relatively high, reaching approximately 1 km. As shown in Fig. 6a, from 25–27 December, cloud cover was minimal, yet the total downward radiation flux decreased daily, which may be associated with the radiative forcing of aerosols. At noon on 24 December, aerosols were primarily concentrated near the surface, potentially favoring the onset of the stove effect. However, the radiosonde temperature profile on the morning of 24 December (Fig. 6f) indicates the presence of a near-surface inversion, which may delay surface warming and the development of turbulence in the lower atmosphere and could partly explain the relatively small variation in surface PM_2.5 concentrations. Starting at 12:00 CST (China Standard Time, UTC+8) on 25 December, a pronounced aerosol layering was observed around 1.1 km. This layer persisted for 12 h with a thickness of 200 m. After sunset, surface radiative cooling led to a decrease in near-surface saturation vapor pressure and an increase in relative humidity, providing favorable conditions for aerosol hygroscopic growth and accelerating liquid-phase and heterogeneous reactions (Cheng et al., 2016; Tie et al., 2017). Consequently, both surface PM_2.5 concentrations and aerosol optical parameters increased significantly during the early hours of 26 December.

From noon on 28 December, cloud cover reduced downward radiation fluxes, allowing aerosols to continue accumulating near the surface. During this cloud episode, the occurrence of virga (in Fig. 6a and c) further enhanced aerosol hygroscopic growth, which could contribute to haze formation and the intensification of pollution. In the afternoon of 29 December, a brief urban heating effect was observed, which could partially enhance turbulent diffusion. However, the limited duration of solar radiation and relatively weak horizontal advection restricted the removal of aerosols. At night, surface pressure decreased, leading to pronounced convective subsidence. Under conditions of high relative humidity, PM_2.5 concentrations and aerosol optical parameters reached their peak values.

Furthermore, three persistent “capping inversion” layers with a dome-like structure were observed during the pollution development stage and are hereafter referred to as dome-type TIs, as indicated by the green dashed regions in Fig. 6d. The identification of dome-type TIs relies primarily on the geometric structure of the vertical temperature profiles and the corresponding aerosol-layer stratification. A comparison with Fig. 8o shows that the daily mean surface PM_2.5 concentrations tend to increase as the dome height decreases and decrease as the dome height increases. Although the dome height on 28 December was relatively high, severe pollution was still observed, suggesting that factors such as aerosol hygroscopic growth under high relative humidity may have contributed.

The hourly evolution of the aerosol vertical structure, atmospheric temperature, and buoyancy acceleration on 25 December is shown in Fig. 10a. Figure 10b presents the variations in temperature tendency, temperature gradient, and TI depth of the elevated TI layer. Figure 10c shows the aerosol heating rate, which is approximately estimated using a simplified shortwave heating formulation in the vertical direction ($Q=\mathrm{1}/\left(\mathit{\rho }{c}_{\mathrm{p}}\right)\cdot \mathrm{d}F/\mathrm{d}z$), where ρ and c_p denote the air density and the specific heat capacity of air, respectively, and F represents the net radiative flux. The lidar ratio at 355 nm ($S=\mathit{\alpha }/\mathit{\beta }$) is shown in Fig. 10d.

Figure 10 Evolution of the vertical structure and atmospheric buoyancy on 25 December. Panels (a)–(d) show the aerosol and temperature structure, buoyancy acceleration, aerosol heating rate, and lidar ratio (355 nm), respectively. In panel (a), the black line with circle markers represents the aerosol backscatter coefficient, and the blue dashed line denotes the temperature profile. Shading indicates buoyancy acceleration, with blue (red) corresponding to negative (positive) values and downward (upward) buoyant motion. Upward- and downward-pointing triangles indicate the TI layer. In panel (b), the green dotted line shows the temperature gradient, while the blue dash–dotted line with square markers represents the TI depth.

Under the combined influence of the elevated TI and weak surface convection, a dome-like stratified structure began to develop near the UABL. As illustrated in Fig. 10b, after 10:00 CST, the temperature tendency around 1.2 km increases significantly and is markedly higher than that in the lower layers, whereas the aerosol heating rate in Fig. 10c remains relatively weak. Backward air-mass trajectory analysis shows that the airflow at approximately 1.5 km is advected from southern regions toward the observation site, indicating a northward transport of warm air. Consequently, horizontal advection may play an important role in the pronounced warming near the upper boundary layer (UABL). Under the combined influence of advective warming and aerosol radiative heating (in Fig. 10c and d), the elevated TI appears to be further intensified. This process is likely to suppress the upward transport of turbulence and induce local aerosol subsidence. Meanwhile, surface warming may enhance buoyancy in the lower atmosphere (in Fig. 10c and d), promoting turbulent uplift of near-surface aerosols. However, the heating rate near the surface remains relatively weak.

As a pronounced stratification develops near the UABL, the heating rate within the aerosol-stratified layer increases substantially, while the warming rate in the lower layer tends to decrease. This reduction may be partly attributed to the attenuation of downward radiative fluxes by the stratified aerosol layer, which could constrain surface warming. Given the mid-latitude location of the study region, the limited duration of daytime solar radiation may have been insufficient to sustain continuous surface heating. Consequently, subsidence within the upper aerosol layer likely became the dominant mechanism regulating vertical aerosol mixing during this period. These interpretations are based on observational consistency and are subject to uncertainties associated with the simplified heating rate estimation.

As shown in Fig. 10b, the temperature gradient and TI depth exhibited nearly opposite variations. Between 12:00 and 15:00 CST, the temperature gradient within the elevated TI increased rapidly from approximately 15 K km⁻¹ to about 37 K km⁻¹, while the TI depth decreased from around 550 m to roughly 240 m. This phenomenon was likely associated with the dome effect induced by aerosol stratification. After 16:00 CST, the temperature gradient decreased markedly, accompanied by a corresponding increase in TI depth, which may have resulted from enhanced vertical mixing of aerosols. During nighttime, variations in both temperature gradient and TI depth were primarily governed by surface longwave radiation and the radiative forcing of aerosols. Therefore, the dome effect played a pivotal role during the aerosol accumulation stage. The negative buoyancy acceleration induced by the TI was the dominant factor suppressing vertical diffusion of aerosols, causing the convective mixing height during the mid-pollution stage to remain confined below approximately 1 km. Although direct measurements of aerosol absorption (e.g., from an aethalometer) were not available during this experiment, the lidar ratio can provide indirect information on aerosol optical properties. As shown in Fig. 10d, the increased lidar ratio in the lower layer during the pollution period may suggest enhanced aerosol absorption and support the development of the near-surface stove effect, although this inference has not been directly verified by dedicated absorption measurements.

3.2.2 Pollution dissipation stage

The lidar and radiosonde temperature profiles indicate that a surface-based TI layer with a thickness of up to 600 m developed in the morning of 30 December, effectively suppressing the vertical diffusion of aerosols. Consequently, despite strong north-westerly winds exceeding 8 m s⁻¹ above approximately 1 km during the daytime (Fig. 6h), surface PM_2.5 concentrations remained persistently high (Fig. 4a). As surface heating intensified, the associated decrease in relative humidity may have promoted aerosol drying and enhanced turbulent mixing in the lower atmosphere, leading to a noticeable reduction in near-surface PM_2.5 concentrations. However, the duration of surface heating was relatively short, potentially limiting the sustained development of turbulence and thereby restricting efficient vertical transport of aerosols. Later in the day, PM_2.5 concentrations increased again, possibly associated with aerosol subsidence and hygroscopic growth.

In the early morning of 31 December, surface PM_2.5 concentrations decreased markedly, likely due to the strong northeasterly winds below 1 km, with wind speeds exceeding 12 m s⁻¹. Although the near-surface PM_2.5 concentrations decreased significantly, this was not attributed to a genuine clean-air advection. Instead, the reduction occurred because nighttime emissions were relatively weak (Le et al., 2020), causing the northeasterly winds–typically associated with polluted air masses–to become relatively clean, thereby diluting the surface PM_2.5 concentrations. At around 14:00 CST on 31 December, a localized stratification was observed, with a formation mechanism similar to that on 25 December. The pollution dissipation during this episode was incomplete. The dissipation on 2 January 2024 followed a similar mechanism to the first event, but was associated with stronger northwesterly winds. These winds not only enhanced horizontal dispersion but also promoted vertical mixing of aerosols. After noon, as near-surface aerosols underwent drying, rapid surface warming increased buoyancy, facilitating upward transport of aerosols. Meanwhile, strong northwesterly winds exceeding 10 m s⁻¹ above approximately 0.5 km efficiently advected the uplifted aerosols downstream. Given that the upwind northwestern region is sparsely populated with low emission intensity and relatively clean air, such winds are often regarded as “clean winds”. Consequently, the combined effects of radiative heating and strong northwesterly flow resulted in a more complete pollution dissipation during this event.

Figure 11a and b show the hourly evolution of aerosol vertical structure, atmospheric temperature, and buoyancy acceleration, while Fig. 11c and d present the heating rate. The top height of the surface-based TI decreased from approximately 0.7 km at 10:00 CST to about 0.4 km at 13:00 CST. Meanwhile, the near-surface aerosol backscatter coefficient showed an overall decreasing trend. A pronounced gradient in the backscatter profile corresponded well with the inversion top height. Considering the relatively high near-surface relative humidity (up to 80 %) during this period, the observed decrease in backscatter is likely dominated by aerosol drying in the lower layer, possibly accompanied by moderate aerosol radiative forcing. By 13:00 CST, near-surface aerosols were mainly concentrated below 0.3 km, promoting the development of a stove effect. The subsequent vertical aerosol distribution followed this pattern. From 14:00–17:00 CST, aerosols below 1 km exhibited well-developed vertical mixing, as indicated by the temperature and backscatter coefficient profiles in Fig. 11a. The vertical evolution before 12:00 CST on 2 January 2024 resembled that on 30 December, with near-surface aerosols gradually drying and showing a downward tendency. When near-surface aerosols accumulated, the lower-layer heating rate increased markedly, accelerating the dissipation of the surface-based TI (Fig. 11c at 16:00 CST). As depicted in Fig. 11b, the near-surface backscatter coefficient decreased at 16:00 CST, while the upper-layer coefficient increased, indicating substantial vertical mixing of aerosols. Furthermore, the potential temperature profiles during the afternoons of 30 December and 2 January showed little variation with height, and in some cases decreased with height, implying an unstable atmosphere conducive to vertical convection. Vertical velocity measurements indicate that during the early and mid-stages of pollution, near-surface vertical velocities were relatively small and dominated by subsidence, whereas upward motion prevailed during the late stage of pollution.

Figure 11 Hourly evolution during the pollution dissipation stage. Panels (a) and (b) show the vertical structures of aerosols and temperature, respectively, while panels (c) and (d) depict buoyancy acceleration. In panels (a) and (b), the black circular line denotes the backscatter coefficient, while the blue dashed line and the shaded areas represent the temperature profile and buoyancy acceleration, respectively. The blue (red) shading indicates negative (positive) buoyancy acceleration, corresponding to downward (upward) buoyant motion. Upward- and downward-pointing triangles are used to indicate the TI layer.

3.2.3 Statistical Characteristics

The observed pollution episode was governed by the combined influence of multiple factors. During the pollution development stage, TIs constrained vertical convection, acting as the primary barrier to aerosol vertical diffusion. Cloud radiative forcing further weakened lower-atmosphere turbulence, thereby constraining convective uplift. Aerosols, in turn, modified the atmospheric thermal structure and buoyancy through radiative interactions involving solar shortwave and surface longwave radiation, leading to complex feedbacks with meteorological conditions. During the pollution dissipation stage, sufficient solar radiation enhanced atmospheric buoyancy, facilitating vertical aerosol mixing, while efficient horizontal advection supported aerosol removal.

Hygroscopic growth further amplified particle size and optical effects, intensifying pollution. Statistical analysis combining near-surface relative humidity and PM_2.5 concentrations (in Fig. 12) reveals a strong positive correlation below 200 m, indicating that high humidity exacerbates aerosol accumulation. Additionally, the vertical gradient of relative humidity below 300 m decreases with increasing PM_2.5, likely due to surface-based TIs and weak turbulence that limit vertical moisture transport during pollution episodes.

Figure 12 Relationships between PM_2.5 concentrations and (a) relative humidity, and (b) its vertical gradient.

Figure 13a presents the vertical distributions of the UABL, the elevated TI, and the surface-based TI. Figure 13b and c show the mean vertical temperature gradients of the elevated and surface-based TIs along with the corresponding PM_2.5 concentrations. It is evident that the elevated TI is closely associated with the UABL. Both TIs exhibit distinct diurnal variations in the mean temperature gradient, with the mean gradient of the elevated TI negatively correlated with PM_2.5 concentrations, whereas that of the surface-based TI shows a positive correlation. Figure 14 presents the frequency distribution of the mean temperature gradient of the elevated TI and the correlation between the TI-top height and the UABL. The results show that TIs with an average temperature gradient ranging from 5 to 10 K km⁻¹ occur most frequently (36 %), and the TI-top height exhibits a significant positive correlation with the UABL height (R²≈0.57).

Figure 13 TI characteristics and aerosol boundary layer distribution. (a) Vertical distribution of the UABL and TIs. The black solid line indicates the UABL. The red shaded area represents the elevated TI layer, with red upward-pointing triangles and blue downward-pointing triangles indicating the top and base heights of the elevated TI, respectively. The blue shaded area represents the surface-based TI layer, with orange upward-pointing triangles and black downward-pointing triangles indicating the top and base heights of the surface-based TI, respectively. (b, c) Black dashed lines represent PM_2.5 concentrations, while red and blue solid lines show the mean vertical temperature gradients of the elevated and surface-based TIs, respectively.

Figure 14 Characteristics of the elevated TI. (a) Frequency of temperature gradient occurrence; (b) Relationship between the UABL and the top height of the elevated TI.

The mean temperature gradient of the elevated TI exhibits a piecewise negative correlation with PM_2.5 concentrations (in Fig. 15a). During intensified pollution, aerosols were primarily concentrated near the surface, enhancing daytime radiative heating in the lower atmosphere and thereby weakening the elevated TI. The observed piecewise pattern may be associated with hygroscopic growth of aerosols during the mid-pollution stage. In contrast, the mean temperature gradient of the surface-based TI shows a pronounced positive correlation with PM_2.5 concentrations (in Fig. 15b).

Figure 15 Relationships between PM_2.5 concentrations and the temperature gradients of (a) the elevated TI and (b) the surface-based TI.

Figure 16 Diurnal rates of change for (a) PM_2.5 concentrations, (b) elevated TI temperature gradient, and (c) surface-based TI temperature gradient.

The diurnal rates of change for PM_2.5 concentrations, the temperature gradient of the elevated and surface-based TIs are shown in Fig. 16. The diurnal rate of change represents the temporal rate of variation within a day, where positive values denote an increase and negative values denote a decrease. PM_2.5 concentrations showed a slight decrease from 00:00 to 10:00 CST, followed by a marked decline between 12:00 and 15:00 CST, and then increased rapidly after 17:00 CST. The elevated TI typically intensified in the early morning, with the strongest rate of change in temperature gradient occurring around 08:00 CST, likely associated with residual aerosol layers from the preceding night. Around noon, enhanced solar radiation warmed both the surface and the top of the aerosol layer, resulting in relatively minor variations in the elevated TI gradient. As solar radiation persisted, enhanced vertical mixing of aerosols contributed to the weakening or disappearance of the elevated TI. Surface-based TI generally strengthened during the early morning and after sunset. The rate of change in surface PM_2.5 was strongly coupled with that of the surface-based TI, with stronger TIs corresponding to higher PM_2.5 concentrations and weaker TIs to lower ones. Notably, TI peaks tended to precede PM_2.5 peaks, reflecting the time lag required for aerosol accumulation. Two distinct minima in PM_2.5 concentrations were observed–one around 04:00 CST and another around 14:00 CST. The nighttime minimum was primarily due to reduced anthropogenic emissions, while the midday minimum occurred when solar radiation was strongest and both surface-based and elevated TI weakened, enhancing vertical mixing and resulting in the lowest near-surface aerosol concentrations. In contrast, PM_2.5 concentrations peaked from late afternoon to nighttime, primarily due to two factors: the strengthening of surface-based TIs and the reduction in solar radiation. These processes suppressed vertical mixing and facilitated the rapid accumulation of aerosols near the surface.

Figure 17 Aerosol accumulation and dissipation during the pollution development. The gray shading indicates the aerosol layer. Orange downward arrows represent downward radiative flux, and green upward arrows indicate the reflected component. The red solid line denotes the temperature profile.

4 Conclusions

Using a self-developed Raman–Mie scattering lidar system from Xi'an University of Technology, the vertical structure of the boundary layer during a severe winter haze event in Xi'an was investigated. A calibrated retrieval algorithm was applied to derive the profiles of temperature, relative humidity, and aerosol, providing high-resolution spatiotemporal observations of the vertical structure during the haze episode. By integrating collocated radiosonde and surface meteorological data, the key meteorological characteristics, influencing factors, and interaction mechanisms governing the formation and evolution of this haze event were analyzed. Lidar observations revealed complex vertical distributions of temperature, relative humidity, and aerosols throughout the pollution episode. The results clearly demonstrated the presence of the aerosol dome effect and hygroscopic growth, both contributing to the intensification of pollution, as well as the mechanisms responsible for its dissipation. The coexistence of elevated and surface-based TIs was a dominant meteorological feature during the haze event, with TI layers and aerosol stratification occurring simultaneously, forming a multilayered aerosol structure.

During the haze development stage, at least three dome phenomena and three cloud-layer episodes were identified. These features weakened vertical mixing within the boundary layer, thereby promoting and intensifying aerosol accumulation. The surface-based TI exhibited a pronounced diurnal variation, serving as another key factor constraining haze dispersion. Notably, the height of UABL did not decrease significantly, with a minimum value of 0.9 km. Within the surface-based TI, aerosols displayed a relatively uniform vertical distribution, suggesting that the TI layer effectively limited the variability of the height of UABL. During daytime, radiative forcing associated with different aerosol vertical distributions primarily manifested as upper-layer or lower-layer heating, altering the vertical temperature gradient. In the absence of strong dynamical processes, rapid aerosol subsidence at night could enhance near-surface warming during the following day, facilitating aerosol vertical mixing. However, the nighttime evolution of the aerosol vertical structure was strongly constrained by surface radiative cooling, which was further modulated by cloud cover and other meteorological factors, rendering the process highly uncertain and complex.

Figure 17 illustrates the main aerosol accumulation and dissipation processes during this pollution episode. Under calm and stable synoptic conditions, for a decreasing-type vertical aerosol profile, when the height of UABL exceeded 1.2 km, aerosol radiative forcing intensified the elevated TI at the top of the aerosol layer, delaying surface warming. The vertical extent of the TI directly influenced atmospheric buoyancy acceleration, resulting in a dome-like stratified structure near the top of the aerosol layer (1 km), with a thickness of 200 m and a duration of up to 12 h. As surface radiative heating progressed, buoyancy in the lower atmosphere increased, driving turbulent uplift of near-surface aerosols and promoting vertical mixing. However, the extent of this mixing was limited by the duration of surface heating. Subsidence of aerosols induced by nighttime radiative cooling was the primary mechanism breaking the stratification. When aerosols were concentrated near the surface, aerosol radiative heating further accelerated near-surface warming, enhancing turbulent diffusion. Statistical analysis revealed that PM_2.5 concentrations were positively correlated with both near-surface relative humidity and the vertical temperature gradient of the surface-based TI. Furthermore, the UABL height showed a strong positive correlation with the top height of the TI layer (R²≈0.57). For the TI layer near the UABL, the temperature gradient exhibited a piecewise negative correlation with PM_2.5 concentrations. Diurnal variations were also evident for PM_2.5 concentrations as well as the temperature gradients of elevated and surface-based TIs. Specifically, the temperature gradient of the elevated TI increased during daytime and decreased at night, consistent with radiative heating at the top of the aerosol layer. In contrast, the temperature gradient of the surface-based TI decreased after sunrise and increased after sunset, reflecting the influence of surface radiative forcing.

This study employed lidar observations to investigate the vertical structure of a severe haze event. However, as the measurements were limited to a single site, it should be noted that haze episodes typically occur over a broader region. The vertical structure of haze may vary due to factors such as local topography and synoptic conditions, potentially giving rise to different phenomena. Comprehensive, large-scale observational networks are therefore required to capture these spatial variations, and such studies will be pursued in future research.