Sensitivity Analysis of the Differential Atmospheric Transmission in Water Vapour Mixing Ratio Retrieval from Raman Lidar Measurements

1. Introduction

Water vapour is a crucial and highly variable greenhouse gas in the Earth’s atmosphere [1], which can significantly influence radiative balance, energy transport [2,3,4], and photochemical processes [5]. It can also affect the radiative budget indirectly through cloud formation and by altering the size, shape, and chemical composition of aerosol particles [6,7,8]. Therefore, systematic and accurate observations of water vapour are essential to understand its impacts and improve climate projections, particularly with high temporal and vertical resolution in the lower atmosphere [9,10,11,12,13,14,15]. Radiosondes (RS) have been used as a benchmark for high vertical resolution water vapour profiles for decades. However, despite substantial advancements in radiosonde technology, they are affected by many uncertainties such as the time to ascend (∼1h) or the dry bias during daytime [16], which must be characterised and corrected. Furthermore, RS launches are spatially dispersed and have low temporal resolution, which depends on the launch frequency, which is often insufficient for many meteorological applications [16]. In this sense, Raman lidar has demonstrated its ability to capture the spatial and temporal evolution of water vapour in the lower atmosphere [2,9,10,11,12,13,14,15].

Raman lidar systems retrieve the water vapour mixing ratio (r) from the intensity of Raman scattering by water vapour ($H_{2}O$) relative to a molecular reference (MR) signal, which can be provided by nitrogen ($N_2$), oxygen ($O_2$), or a combination of both molecules. Accurate retrievals require correction for temperature–dependent Raman cross–sections, a system–specific calibration constant (K), and a differential atmospheric transmission $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation, which accounts for extinction differences between the Raman–shifted wavelengths of molecular reference (${\lambda }_{MR}$) and water vapour (${\lambda }_{H_{2}O}$). While the temperature dependence was analysed in detail by [17], with the required filter passbands determined to minimise temperature sensitivity in the troposphere, the effect of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ has not yet been thoroughly investigated. Differences in $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ arise primarily from Rayleigh scattering by air molecules, which depends on ${\lambda }^{-4}$ and can be calculated from temperature and pressure profiles. However, calculations of the aerosol contribution are more complex due to their large spatial and temporal variability and the general difficulty of measuring aerosol extinction near the ground using lidar. The first lidar systems for water vapour measurements used the Raman vibrational–rotational (VR) configuration with $N_2$ and $H_{2}O$ channels centred at ∼387 nm and ∼408 nm. The transmission ratio between Raman–shifted wavelengths of nitrogen and water vapour can deviate by 5-10% from unity, and depending on aerosol load, may introduce an additional 5% [2,18]. However, the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation is still often neglected in many studies e.g., [11,13,19,20,21,22,23].

More recently, some Raman lidar systems–such as those integrated into the European Aerosol Research Lidar Network (EARLINET), which is part of the Aerosol, Clouds, and Trace Gases Research Infrastructure ACTRIS, [24], and the Latin American Lidar Network LALINET, [25]–which have traditionally focused on aerosol studies—have implemented detection of pure–rotational Raman (RR) signals. These signals provide significantly higher Signal–to–Noise Ratio (SNR), enabling shorter integration times and even daytime Raman retrievals [26,27] and atmospheric temperature retrievals [28,29]. In such cases, neglecting the differential atmospheric transmission term may lead to larger uncertainties in the retrieval of the calibration constant for water vapour, and thus in the water vapour mixing ratio measurements. In this configuration, both molecular and aerosol scattering may become important due to the significant wavelength difference between molecular reference, centred at ∼354 nm, and VR water vapour, centred at ∼408 nm, with a spectral difference of ∼54 nm.

Not accounting for the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ effects introduces systematic uncertainties in water vapour retrievals. However, if these effects are appropriately characterised and calculated, such uncertainties can be transformed into random ones. Actually, random uncertainties can be treated statistically and reduced by increasing the number of observations, as their expected value is zero by definition [30,31,32]. Therefore, a proper characterisation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ can enhance the measurements of Raman lidars.

This study presents a sensitivity analysis of the differential atmospheric transmission term $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ in Raman lidar water vapour retrievals. The sensitivity study focusses on two different configurations for measuring water vapour mixing ratio using Raman lidar measurements, the vibrational–rotational nitrogen at ∼387 nm and the pure–rotational signal at ∼354 nm for the molecular reference, both with water vapour at ∼408 nm and thus the latter presenting larger spectral differences when compared with the former. The impact of molecular and aerosol extinction on the differential transmission term is investigated in detail. In particular, aerosol extinction is retrieved from Raman lidar measurements performed at the Granada urban station (UGR, 37.16ºN, 3.60ºW, 680 m a.s.l.) within the framework of EARLINET/ACTRIS, and it is also simulated using an exponential decay approach based on Aerosol Optical Depth (AOD) from co–located Aerosol Robotic Network AERONET, [33] sun photometer (SP) observations. The AOD–based approach enables the operational calculation of differential transmission, avoiding the considerable effort required to retrieve aerosol extinction profiles operationally using the lidar technique, allowing the analysis of long–term Raman lidar water vapour datasets. A statistical comparison is carried out between the simulated atmospheric transmission calculation and those retrieved by the Raman lidar.

The structure of the paper is as follows: Section 2 describes the experimental station and the instrumentation used in the study. Section 3 provides an overview of the Raman lidar technique for obtaining water vapour measurements, along with a detailed explanation of the approach used to calculate the differential atmospheric transmission. Section 4 presents a discussion of the main results, while Section 5 outlines the key conclusions and suggestions for future work.

2. Experimental Site and Instrumentation

2.1. Granada Urban Station

The experimental part of the study was recorded at the Andalusian Global ObseRvatory of the Atmosphere (AGORA) of the University of Granada. AGORA is located in Granada (37.16ºN, 3.60ºW, 680 m above sea level (asl), a medium–size lightly–industrialized city situated in a natural basin surrounded by mountains with altitudes over 1000 m. The climate exhibits continental–Mediterranean characteristics, with significant seasonal temperature variations–cool winters and hot summers. The geographical location, being relatively close to the African continent (about 200∼km) and approximately 50∼km from the western Mediterranean basin, also implies that the station is frequently influenced by Saharan dust intrusions [34].

2.2. Instrumentation

Lidar measurements were performed using the Raman lidar system MULHACEN (Figure 1a) based on a LR331D400 manufactured by Raymetrics S.A., Greece. This system is located at the Granada urban station and has been part of the EARLINET network since April 2005 [11] and contributes to the ACTRIS research infrastructure [24]. The lidar system is configured in a monostatic biaxial alignment pointing vertically to the zenith. The light source is a Nd:YAG pulsed laser (Quantel CFR Series), with a fundamental wavelength at 1064 nm (110 mJ), and additional emissions at 532 nm (65 mJ) and 355 nm (60 mJ) using harmonic generators, with a repetition rate of 10 Hz and pulse duration of 8 ns. The radiation is collected by a 400 mm diameter Cassegrain telescope. The receiving subsystem also includes a wavelength separation unit with dichroic mirrors, interference filters, and a polarization cube. The instrument operates with a vertical resolution of 7.5 m. Due to the instrument setup, incomplete overlap between the laser beam and the receiver field of view limits the observations near the surface [35]. Detection occurs in seven channels corresponding to elastic wavelengths at 1064 nm, 532 nm (parallel and perpendicular polarized), and 355 nm, as well as inelastic wavelengths at ∼607 nm (nitrogen Raman–shifted signal excited by radiation at 532 nm), ∼387 nm (nitrogen, Figure 1b) and ∼408 nm (water vapour) Raman–shifted signals, both excited by 355 nm radiation. In May 2017, a new configuration of dichroic mirrors based on rotational Raman filters were implemented to ensure sufficient signal–to–noise ratio for aerosol retrievals [26], and as a result, the pure rotational Raman channel for the molecular reference at ∼354 nm replaced the vibrational–rotational nitrogen channel at ∼387 nm [27]. This change affected the optical path for the detection of both molecular reference and water vapour (Figure 1c), impacting the overlap functions, the differential atmospheric transmission calculations and, consequently, the calibration constants used to retrieve the water vapour mixing ratio. All modifications were carried out to optimize the system performance and comply with the EARLINET/ACTRIS quality standards for aerosol retrievals [36].

The AGORA station also operates a SP within the framework of the AERONET [33]. The main instrument of AERONET is the Cimel sun photometer, which performs direct sun and sky radiance measurements at various wavelengths and provides AOD at 340, 380, 440, 670, 870, and 1020 nm. Particularly, we used Level 2.0 version 3 data [37]. In addition, the AGORA station is equipped with a GRAW DFM–09 RS, a lightweight radiosonde that provides measurements of temperature (resolution 0.01 ºC, accuracy 0.2 ºC), pressure, (resolution 0.1 hPa, accuracy 0.3 hPa), relative humidity (RH, resolution 1%, accuracy better than 4%) [13].

3. Methodology

3.1. Water Vapour Mixing Ratio Retrieval from Raman Lidar Measurements

The water vapour Raman lidar technique is based on the ratio of Raman scattering intensities from water vapour and the molecular reference [2,9]. The Raman lidar equation for the molecular reference and water vapour signals can be written as follows:

$$P(z,{\lambda }_i)=P\left({\lambda }_0\right)·K_i·{\displaystyle \frac{O_i\left(z\right)}{z^2}}·F_i\left[T\left(z\right)\right]N_i\left(z\right)[d{\sigma }_i({\lambda }_0,\pi )/d\Omega ]·exp\left\{-{\int }_0^z[\alpha (z^{\prime },{\lambda }_0)+\alpha (z^{\prime },{\lambda }_i)]d{z}^{\prime }\right\}$$

(1)

where index i indicates the molecular reference or water vapour molecule; $P(z,{\lambda }_i)$ is the backscattered signal from range z at the Raman–shifted wavelength ${\lambda }_i$; $P\left({\lambda }_0\right)$ is the emitted laser power at the wavelength ${\lambda }_0$; $K_i$ is the range-independent system calibration constant, which accounts for the system optical efficiency, the telescope receiver area, the photomultiplier tube (PMT) spectral efficiency, and the laser output energy. $O_i$ is the overlap function; $F_i\left[T\left(z\right)\right]$ is the temperature–dependent function of the Raman scattering [9,17]; $N_i\left(z\right)$ is the number density; $d{\sigma }_i({\lambda }_0,\pi )/d\Omega$ is the cross section at the Raman–shifted wavelength; and the exponential term represents the atmospheric transmission function ($\mathcal{T}$) due to the total extinction coefficient $\alpha$ at wavelengths ${\lambda }_0$ and ${\lambda }_i$.

The ratio $N_{H_{2}O}\left(z\right)/N_{MR}\left(z\right)$ is proportional to the water vapour mixing ratio (r), which is defined as the mass of water vapour divided by the mass of dry air in a sample of the atmosphere [39] allowing the Raman lidar technique to provide a direct measurement of the atmospheric water vapour [2,9]. Assuming identical overlap functions and range–independent Raman Backscatter cross sections for the molecular reference and $H_{2}O$:

$${\displaystyle \frac{N_{H_{2}O}\left(z\right)}{N_{MR}\left(z\right)}}={\displaystyle \frac{P(z,{\lambda }_{H_{2}O})}{P(z,{\lambda }_{MR})}}·{\displaystyle \frac{F_{MR}\left[T\left(z\right)\right]}{F_{H_{2}O}\left[T\left(z\right)\right]}}·{\displaystyle \frac{K_{MR}{\sigma }_{MR}}{K_{H_{2}O}{\sigma }_{H_{2}O}}}·exp\left\{{\int }_0^z[\alpha (z^{\prime },{\lambda }_{H_{2}O})-\alpha (z^{\prime },{\lambda }_{MR})]d{z}^{\prime }\right\}$$

(2)

and thus,

$$r\left(z\right)={\displaystyle \frac{P(z,{\lambda }_{H_{2}O})}{P(z,{\lambda }_{MR})}}·{\displaystyle \frac{F_{MR}\left[T\left(z\right)\right]}{F_{H_{2}O}\left[T\left(z\right)\right]}}·K·exp\left\{{\int }_0^z[\alpha (z^{\prime },{\lambda }_{H_{2}O})-\alpha (z^{\prime },{\lambda }_{MR})]d{z}^{\prime }\right\}$$

(3)

where the term $P(z,{\lambda }_{H_{2}O})/P(z,{\lambda }_{MR})$ represents the backscattered signal ratio, $F_{MR}\left[T\left(z\right)\right]/F_{H_{2}O}\left[T\left(z\right)\right]$ is the ratio of the temperature–dependent functions for the Raman molecular reference and $H_{2}O$ channels, K is the calibration constant. The exponential term represents the differential atmospheric transmission $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ that occurs at the two Raman–shifted wavelengths, ${\lambda }_{MR}$ and ${\lambda }_{H_{2}O}$, where the total extinction coefficient $\alpha (z^{\prime },\lambda )$ at range z and wavelength $\lambda$ is the sum of the molecular extinction ${\alpha }_{mol}(z^{\prime },\lambda )$ and aerosol extinction ${\alpha }_{aer}(z^{\prime },\lambda )$ coefficients. The equation can be rewritten as:

$$r\left(z\right)={\displaystyle \frac{P(z,{\lambda }_{H_{2}O})}{P(z,{\lambda }_{MR})}}·{\displaystyle \frac{F_{MR}\left[T\left(z\right)\right]}{F_{H_{2}O}\left[T\left(z\right)\right]}}·K·\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$$

(4)

The assumption of identical overlap functions for the molecular reference and $H_{2}O$ might not be true in real applications and differences between both overlap functions are found in the near range [9]. The temperature dependence of pure rotational and vibrational–rotational Raman was presented by [17], who assessed the required filter passbands for water vapour measurements with minimal temperature sensitivity using Nd:YAG lasers. According to this assessment the filters used in the MULHACEN Raman lidar for water vapour (centred at ∼408 nm, Full Width at Half Maximum (FWHM) = 1 nm) and molecular reference (centred at ∼387 nm, FWHM = 2.7 nm, or centred at ∼354 nm, FWHM = 0.8 nm) select the Raman Stokes spectrum (∼387 nm) or Raman Anti–Stokes spectrum (∼354 nm), resulting in essentially temperature independent measurements of the Raman signals, and making it reasonable to approximate the ratio $F_{MR}\left[T\left(z\right)\right]/F_{H_{2}O}\left[T\left(z\right)\right]$ to 1 [17]. Once the differential atmospheric transmission is accounted for, the calibration constant (K) is determined by comparing lidar profiles with simultaneous and co–located radiosonde measurements of water vapour using the iterative least squares method [13].

3.2. Differential Atmospheric Transmission Assessments

The differential atmospheric transmission $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ in the retrieval of the water vapour mixing ratio arises from the difference in atmospheric transmission at the two wavelengths, ${\lambda }_{MR}$ and ${\lambda }_{H_{2}O}$. This difference is mainly due to the ${\lambda }^{-4}$ dependence of Rayleigh scattering by air molecules [2], with an additional contribution from the wavelength–dependent Mie scattering by aerosols. The molecular extinction coefficient can be calculated for Rayleigh scattering [40,41] using temperature $T\left(z\right)$ and pressure $P\left(z\right)$ profiles. In this study, these profiles were obtained from ECMWF model data for the location of Granada [42], available from the ACTRIS Data Centre: https://hdl.handle.net/21.12132/1.16d392060df54287. The aerosol extinction coefficient can be retrieved directly using the Raman inversion method [43], or indirectly from aerosol backscatter coefficients retrieved with the Klett inversion method [44,45], multiplied by an assumed lidar ratio. In this study, the aerosol extinction coefficient was obtained by multiplying the aerosol backscatter coefficient (${\beta }_{aer}$) at 355 nm by a height–independent lidar ratio of 50 sr. This lidar ratio value has been used in previous studies at our station [46,47,48]. To estimate the aerosol backscatter coefficients at the Raman wavelengths (354, 387, and 408 nm), the spectral dependence of the backscatter was considered using the Ångström law [49], with the Ångström exponent (AE) derived from ${\beta }_{aer}$ at 355 and 532 nm, both of which are emitted by MULHACEN. These profiles are also available from the EARLINET network (https://data.earlinet.org/earlinet/).

However, the operational and automatic retrieval of aerosol backscatter or extinction coefficient requires considerable computational effort. As a result, the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ term is often neglected in water vapour Raman lidar retrievals in several studies [e.g., [11,13,19,20,21,22,23]. The calculation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ is a crucial step in the accurate retrieval of water vapour from Raman lidar measurements, as uncalculated systematic uncertainties on the order of 5–15% can introduce long–term biases in water vapour measurements. In addition, depending on the aerosol load, up to an additional 5% may be required in this term [2,18]. The impact of neglecting this term becomes more significant for Raman lidars that retrieve water vapour using VR water vapour and RR channels, due to the larger spectral separation between the molecular reference and water vapour Raman wavelengths. In contrast, once properly accounted for, these uncertainties can be treated as random. Since random uncertainties have an expected value of zero by definition, their impact can be reduced by increasing the number of observations [30,31,32].

To address this challenge, the present study proposes an approximation for the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation aiming to overcome the considerable computational effort required for the automatic and operational retrieval of aerosol extinction coefficients. This approach is based on AOD values obtained from a sun photometer and the modelling of the vertical distribution of aerosol extinction using an exponential decay function of aerosol extinction with altitude [50]. Accordingly, the vertical profile of aerosol extinction coefficient is expressed as:

$${\alpha }_{aer}(z,\lambda )={\alpha }_{aer}(0,\lambda )·exp\left(-{\displaystyle \frac{z}{H_p}}\right)$$

(5)

where ${\alpha }_{aer}(0,\lambda )$ is the surface extinction and $H_p$ is the scale height. In this study, $H_p$ is taken as the ABL height ($H_{ABL}$), which was determined for Granada urban station by [51].

When using AOD($\lambda$) for the total aerosol load, Equation 5 is given as:

$${\alpha }_{aer}(z,\lambda )={\displaystyle \frac{AOD\left(\lambda \right)·exp\left(-\frac{z}{H_{ABL}}\right)}{{\int }_0^{\infty }exp\left(-\frac{z^{\prime }}{H_{ABL}}\right)d{z}^{\prime }}}$$

(6)

AOD values at 354, 387, and 408 nm were derived from the AOD at 380 nm using the AE between 340 and 440 nm from sun photometer measurements.

4. Results and Discussion

4.1. Sensitivity of the Differential Atmospheric Transmission with Wavelength

The sensitivity of the differential atmospheric transmission to wavelength is shown in Figure 2. The upper panels correspond to the first optical configuration used to retrieve the water vapour mixing ratio with the MULHACEN Raman lidar system, which employed the VR nitrogen channel at ∼387 nm, while the lower panels present results for the second configuration, which used the RR signal at ∼354 nm for the molecular reference. In both configurations, the water vapour Raman signal was measured at ∼408 nm. Particularly, Figure 2 presents the mean and standard deviation (SD) profiles of ${\alpha }_{mol}$ (panels a, d), ${\alpha }_{aer}$ (panels b, e) and $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ (panels c, f) for the period 2010–2016, during which a total of 300 lidar aerosol inversions were performed. Molecular extinction profiles were calculated using temperature and pressure data from the ECMWF model for the Granada urban station. Aerosol extinction profiles were derived by multiplying the aerosol backscatter coefficients measured by the MULHACEN lidar system by a lidar ratio of 50 sr. Panel (c) includes a zoomed inset to better illustrate the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ profiles of the first configuration, allowing for a clearer comparison with the second (panel f). A statistical analysis of the mean calculation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ for different layers is also presented in Table 1.

Results from Figure 2 show an increase with height of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ and of its standard deviation (panels c and f, shaded areas), which is explained by the increased scattering due to molecules and aerosols as the path length increase. For the first configuration (VR nitrogen at ∼387 nm), neglecting the total differential transmission term may introduce systematic uncertainties of approximately 3.4% at 3 km, where most of the atmospheric water vapour is typically concentrated. The uncertainty may increase to 5.1% at 6 km (black curve in Figure 2c). On the other hand, for the second configuration that uses RR signal at ∼354 nm (Figure 2, lower panels) larger values in $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ are observed. These reach approximately 10% at 3km and 15% at 6 km (black curve in Figure 2f) (Table 1). These higher values of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ are explained by the larger spectral difference between the RR signal (∼354 nm) and water vapour (∼408 nm) signal. It should be noted that, the molecular component influence becomes more considerable than the aerosol one due to the ${\lambda }^{-4}$ of molecular scattering as opposed to the roughly ${\lambda }^{-1}$ dependence of aerosol scattering. Specifically, $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ reaches about 7.0% at 3 km and 12% at 6 km (red curve in Figure 2f) . Similarly, the aerosol contribution (blue curve in Figure 2f) presents larger $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ values than in the previous configuration, with values of approximately 3.3% at 3 km and 4.0% at 6 km (Table 1).

Figure 2 and Table 1 clearly show that not including the calculation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ in Equation 3 would introduce a significant systematic effect in water vapour mixing ratio retrievals. According to the recommendations of the Joint Committee for Guides in Metrology in their Guide to the Expression of Uncertainty in Measurement [JCFG/GUM, [30,32], “it is assumed that the result of a measurement has been corrected for all recognized significant systematic effects and that every effort has been made to identify such effects.” This principle is the driving force behind the requirement to include these calculations in the proper evaluation of the water vapour mixing ratio and has been done in previous studies which used the VR nitrogen Raman signal at ∼387 nm [2,18].

We demonstrate that these effects are even more significant when using RR signal at ∼354 nm for the molecular reference. Ideally, ${\alpha }_{mol}$ should be calculated using correlative temperature and pressure profiles from radiosonde measurements. If such data are unavailable, calculations can be made using profiles from numerical weather models, which, despite possible differences, may still be used to convert systematic uncertainties into random uncertainties. For aerosols, the ideal correction uses correlative lidar measurements of ${\alpha }_{aer}$. However, when these are not available, we propose an approach based solely on correlative AOD measurements, which similarly may translate systematic uncertainties into random uncertainties.

4.2. Sensitivity of Differential Atmospheric Transmission on Aerosol Loading

The dependence of the differential atmospheric transmission on aerosol loading is evaluated using different ranges of AOD values, as illustrated in Figure 3, differentiating between the VR (panel a) and RR (panel b) molecular reference configurations. The AOD values were obtained by integrating of the aerosol ${\alpha }_{aer}$ at 355 nm from MULHACEN lidar inversions. Data correspond to the 2010–2016 measurement period at the Granada EARLINET/ACTRIS station, where the mean and SD of AOD₃₅₅ value were of 0.28 ± 0.1. The effect in $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ caused by molecular scattering alone (red curve) is also included for comparison. The aerosol load conditions (blue curves) considered in the calculation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ range from very low AOD values to relatively high aerosol load conditions (solid blue curve). Specifically, low aerosol load corresponds to AOD values below the mean minus one SD (i.e., less than 0.18), medium aerosol load includes AOD values between the mean minus and plus one SD (from 0.18 to less than 0.38, dash–dot blue curve), and high aerosol load refers to AOD values equal to or greater than the mean plus one SD (0.38 or higher). In addition, Table 2 presents a statistical analysis of the mean calculation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$, evaluated across different layers and aerosol loading conditions.

Results from Figure 3 reveal that the influence of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ increases with aerosol loading, showing larger deviations from unity as AOD increases. Consequently, if $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ is neglected, higher aerosol loads lead to greater systematic uncertainties in water vapour mixing ratio retrievals. For the first configuration (panel a), the results show that percentage differences in $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ from unity reach values between 4.5% and 6% up to 6 km for low and high aerosol load conditions, respectively (Table 2). These values are consistent with those reported in previous studies [2,17], and illustrate that $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ must be included in the calculation of mixing ratio in order to convert systematic uncertainties into random ones.

The most relevant result is that the effect of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ is clearly more significant for the second configuration, and this impact strongly depends on aerosol loading. In particular, under locally high aerosol load conditions (solid blue curve in Figure 3b), the total $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation reaches approximately 18% up to 6 km, compared to 12% when considering only the Rayleigh scattering calculation (Table 2). For medium and low aerosol load conditions (dash–dot and dotted blue curves in Figure 3b), the additional mean $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ within the first 6 km required relative to the molecular case are approximately 2% and 3%, respectively (Table 2). Therefore, including $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ in the calculation of water vapour mixing ratio is critical particularly for high aerosol loads that in the UGR station can reach values above 1.0 during extreme Saharan dust events [11,48,52,53].

4.3. Comparison of AOD from Sun Photometer and Raman Lidar Measurements

This study proposes an alternative approach to estimate the aerosol extinction coefficient profiles to overcome the considerable computational effort required for the operational retrieval when the lidar overlap function is significant. The approach is based on SP AOD values and assumes an exponential decay function of aerosol concentration with altitude (Equation 6, Section 3.2). To ensure the quality of the estimated aerosol profiles, SP AOD values were compared with those obtained by integrating aerosol extinction profiles retrieved from MULHACEN Raman lidar measurements. The comparison was carried out using simultaneous daytime AOD values, while for nighttime comparisons, the AOD values were estimated by interpolating between the last value in the previous afternoon and the first value in the next morning, as previously done in other works in the literature e.g., [54].

This evaluation is illustrated in Figure 4 for the period 2010–2016, which presents the total dataset of 300 lidar measurements (panels a, d), daytime measurements (250 cases, panels b, e), and nighttime measurements (50 cases, panels c, f). The upper panels of Figure 4 show a direct inter–comparison of AOD values obtained by lidar versus those obtained by sun photometer, including the 1:1 reference line (dashed black) and the linear regression fit (solid red) with corresponding equations and uncertainties. The colour scales indicate magnitude of the residuals. The lower panels display histograms of the AOD differences (${\mathrm{AOD}}_{Lidar}$-${\mathrm{AOD}}_{SP}$). The parameters of the linear fits, as well as other statistical parameters, are summarized in Table 3. The mean SP AOD values at 355 nm for the total, daytime and nighttime datasets were 0.25 ± 0.09, 0.26 ± 0.09, and 0.23 ± 0.07, respectively.

For the complete and daytime datasets, the linear regressions (Figure 4a, b) show good agreement between lidar and sun photometer ($R^2$ of 0.80 and 0.85, respectively). In both cases, the slope of the linear fits (0.93 ± 0.03 and 0.96 ± 0.03) suggest an overestimation of lidar AOD values compared with reference sun photometer measurements (approximately 7%), which is supported by the bias and root mean square error (RMSE) values of 0.03 ± 0.04 and 0.05, respectively (Table 3). This overestimation may be related to the constant assumption of lidar ratio or to issues associated with the incomplete overlap region. Histograms (Figure 4d, e) display an approximately Gaussian shape centred near zero, although slightly skewed towards negative values, consistent with previous findings indicating a tendency for lidar to overestimate AOD e.g., [55,56,57,58]. Nevertheless, the fact that the histograms are centred close to zero plus the values of the residuals confirm that the majority of differences cluster near zero indicates no significant systematic uncertainties across the AOD range and further supports the robustness of the comparison.

The nighttime comparison (Figure 4c, Table 3) shows somewhat poorer agreement between lidar and sun photometer compared to the previous cases. The linear fit yields an $R^2$ of 0.66, which, although modest, remains statistically significant at the 95% confidence level. The linear regression indicates that for AOD values below 0.2, lidar values tend to be slightly higher than those from the SP, whereas for AOD values above 0.2, lidar values tend to be lower. This pattern is reflected by a small negative bias of −0.01 ± 0.04. The histogram of differences (Figure 4f) is centred at zero but does not show a clear Gaussian shape with a relatively high frequency of values in the extremes. This large dispersion may be attributed to increased uncertainties introduced by the interpolation approach used to estimate nighttime AOD, which could be minimised if actual nocturnal measurements were available, for example, through star [59,60] or moon photometry [61]. Nevertheless, the results demonstrate that estimating nighttime AOD by interpolating daytime sun photometer data provides a reasonable approximation for calculating $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$, because the fundamental requirement for a satisfactory correction is that its mean value be zero [30,32]. If the random component increases through a given procedure, this is acceptable as long as the mean remains zero, and effective can be used to transform systematic uncertainties into random ones albeit with larger random uncertainties during the nighttime.

4.4. Validation of Differential Atmospheric Transmission Estimated from Sun Photometer Against Raman Lidar

The aerosol extinction coefficients profiles are estimated from AERONET sun photometer AOD values using an exponential decay function (Equation 6, Section 3.2). Based on these estimates, $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ due to aerosols is calculated (Eqs. 3 and 4). These simulated $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ are then compared with reference profiles obtained from 300 aerosol extinction profiles retrieved by lidar inversion. Validation of the exponential decay approach is presented in Figure 5, including results for both the first (vibrational–rotational nitrogen at ∼387 nm) and second (pure–rotational signal at ∼354 nm for the molecular reference) optical configurations of the MULHACEN Raman lidar system. Additionally, statistical summaries of the comparisons are provided in Table 4.

For both optical configurations, the simulated $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ profiles–calculated using the exponential decay approach (Figure 5)–closely match the values obtained from the MULHACEN Raman lidar (panels c and g). Within the first 3 km, the mean aerosol contribution to $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ is slightly lower in the simulations, with a 1.1% value compared to 1.2% from lidar inversions for the configuration employing VR nitrogen and water vapour channels (Figure 5c). For the configuration based on RR molecular reference and VR water vapour channels (Figure 5g), the mean value is 3.0% in the simulations and 3.3% in the lidar-derived data. At 6 km, simulated values are slightly lower than lidar-derived ones, with mean values of 1.2% (simulated) compared to 1.3% (lidar), and 3.5% versus 4.0% for the first and second optical configurations, respectively (panels c and g in Figure 5). The results shown in panels d and h of Figure 5 are particularly relevant, as the differences in the aerosol differential transmission calculation between the lidar and simulated values are nearly zero, with relative deviations of 0.10% and 0.40% within the 6 km for the first and second configurations, respectively (Table 4). These findings demonstrate a strong agreement between the simulations and the lidar measurements, with $R^2$ of 1.0, slope of 1.0 ± 0.0, and intercepts of 0.001 ± 0.00 and 0.002 ± 0.00 for molecular reference at ∼387 and ∼354 nm, respectively.

It is important to highlight that the residuals between the lidar-derived and simulated values of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$–using the exponential decay approach–exhibit a mean close to zero and a low standard deviation. This behaviour indicates that, when calculation of differential atmospheric transmission is included that the systematic uncertainties that are introduced when the term $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ is not accounted for are be transformed into random uncertainties. This is particularly relevant for the retrieval of the water vapour mixing ratio, as random uncertainties can be statistically treated and propagated more robustly than systematic ones. Neglecting this term may lead to significant uncertainty in the water vapour retrieval [2,18]. Therefore, the proposed approach can minimise uncertainty and improve the accuracy of water vapour retrievals from Raman lidar measurements. The results demonstrate that this term can and should be calculated automatically and operationally using an exponential decay approach, offering a clear advantage over other studies in which it is often neglected e.g., [11,13,19,20,21,22,23].

Other studies, such as [18], have calculated the term $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ by simulating the aerosol extinction profile using a step function. This approach assumes a well–mixed aerosol layer confined within the Atmospheric Boundary Layer (ABL), characterised by stability and an approximately homogeneous vertical distribution of aerosols. Under these conditions, it is reasonable to assume that the aerosol extinction coefficient remains constant with height within the ABL. Therefore, the total AOD is uniformly distributed throughout the ABL, and the aerosol extinction coefficient is considered constant within this layer and zero above it. We also applied this approach and found that the resulting calculations closely matched those obtained from lidar inversions, with $R^2$ of 1.0, a slope of 1.0 ± 0.0, and intercepts of −0.0001 ± 0.00 and −0.001 ± 0.00 for the molecular reference at ∼387 and ∼354 nm, respectively. The mean relative deviations between the simulated and lidar-derived $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculations within the first 6 km were −0.03% and −0.10% for the first and second optical configurations of the MULHACEN system, respectively. However, the step function approach may be less accurate in the presence of decoupled aerosol layers above the ABL–structures frequently observed at our station [11,48,52,53]. In such scenarios, the exponential decay method offers a more physically realistic representation of vertical aerosol distribution, reflecting natural processes such as the exponential decrease of atmospheric pressure with altitude. Despite their differences, both approaches are useful and provide a well–defined methodology for estimating the profile of aerosol extinction when no correlative measurements are available.

4.5. Evaluation of the Differential Atmospheric Transmission Calculation in Water Vapour Measurements

To evaluate the impact of the differential atmospheric transmission term on the retrieval of water vapour mixing ratio, Raman lidar water vapour profiles were calculated both considering and not considering the term $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ (Equation 3). It was assumed that the temperature dependence of the Raman backscatter cross sections is negligible, since the interference filters used to select the Raman scattering ensure a nearly temperature independent effect [17] (Section 3.1). After calculating (or not calculating) the differential transmission term, the calibration constant K was determined by comparison with radiosonde data, using the robust iterative least squares method proposed by [13]. The calculation $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ (Eqs. 3 and 4) was estimated using an exponential decay function (Equation 6) based on SP AOD data.

Two different cases were analysed to investigate the effect of the differential atmospheric transmission term on Raman lidar water vapour measurements (Figure 6). Case I, on 11 July 2016 at 20:50 UTC (Figure 6, upper panels), corresponds to the period when the MULHACEN Raman lidar used the first optical configuration to retrieve the water vapour mixing ratio (VR nitrogen at ∼387 nm and VR water vapour at ∼408 nm). Case II, on 24 August 2017 at 19:36 UTC (Figure 6, lower panels), represents a study case where water vapour was retrieved using RR molecular reference at ∼354 nm and water vapour at ∼408 nm (second configuration). The Raman lidar water vapour mixing ratio profiles calculated including the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$, using the exponential decay approach are represented by a solid blue curve; and lidar profiles without $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ are indicated by a dotted blue curve. Panels (b) and (e) in Figure 6 display the differences between the two water vapour mixing ratio profiles for Cases I and II, respectively, while the panels (c) and (f) show the relative differences.

Cases I and II correspond to AOD at 355 nm of 0.31 (moderate aerosol load; see Section 4.2) and 0.38 (high aerosol load), respectively. These values result in total differential transmission at 6 km of 4.5% and 12.4%. In both cases, significant differences are observed between the water vapour mixing ratio profiles retrieved with and without calculating the differential transmission term (Figure 6). The maximum relative differences reach approximately −4% in Case I and −11% in Case II (panel c, f). Overall, the profiles excluding the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation exhibit a systematic underestimation of water vapour compared to the ones including the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation, with mean relative differences (using the fully calculated profiles as reference) of −2.3% and −4.0% within the first 6 km (Table 5). The most significant discrepancies are observed in the lower troposphere (below 3 km), where water vapour concentrations are typically highest and there is often strong aerosol concentrations. In this region, mean relative differences of −2.5% and −7.0% are found for Cases I and II, respectively. However, important differences are also present at higher altitudes, where, despite lower absolute humidity, the cumulative effect of differential transmission remains relevant (Table 5).

Furthermore, the impact of the molecular reference channel wavelength selection is clearly evident. In Case II (Figure 6, lower panels), where the spectral separation between the molecular reference and water vapour channels is larger, more pronounced differences are observed compared to Case I (Figure 6, upper panels). For instance, below 4 km, relative differences reach −6.0% in Case II compared to −2.6% in Case I, emphasizing the greater need to account for the differential transmission when the spectral gap increases. Regarding the precipitable water vapour (W) retrieved from lidar, in Case I it is 18.90 ± 0.19 mm with the differential transmission, and 18.42 ± 0.20 mm without it. This corresponds to relative underestimations of 0.7% (with $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$) and 3.3% (without it) compared to the Global Navigation Satellite System (GNSS) value (19.04 ± 0.33 mm), illustrating that the including $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ improves the agreement. In Case II, omitting the calculation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ results in a 6% underestimation, while including it reduces the relative differences to only 0.2% (Table 5, GNSS W = 21.20 ± 0.12 mm). These results are consistent with the mean values reported in Table 1 and Table 2, and highlight the importance of including the differential transmission calculation to prevent systematic uncertainties in both vertical profiles and integrated water vapour estimates. Finally, it is important to note that neglecting this calculation introduces uncertainty in the determination of the calibration constant, as it would be derived from a signal affected by a bias. Consequently, this leads to a systematic uncertainty in the water vapour when using a specific calibration constant for subsequent retrievals. In contrast, our results demonstrate that the proposed approach enables the transformation of this systematic uncertainty into a random one, which can then be statistically treated, ultimately improving the reliability of the final product.

5. Summary and Conclusions

Accurate retrieval of the water vapour mixing ratio (r) from Raman lidar measurements requires precise calculation of the differential atmospheric transmission, $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$, which accounts for the differences in molecular and aerosol extinction between the molecular reference and water vapour Raman–shifted wavelengths. Although this calculation is essential–since neglecting it often leads to significant systematic uncertainties in the retrieval of the water vapour mixing ratio–it has often been neglected in recent studies, potentially introducing long–term biases that compromise climate trend analyses and model validation. These uncertainties become more significant in systems with large spectral separation between Raman channels, such as those using a pure rotational (RR) filter to detect the molecular reference signal at ∼354 nm and a vibrational–rotational (VR) filter at ∼408 nm for water vapour. This is due to the wavelength dependence of molecules and Mie scattering.

This study presents a sensitivity analysis of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ with respect to wavelength and aerosol loading, using lidar inversions and an operationally feasible approach based on sun photometer (SP) aerosol optical depth (AOD) measurements, where aerosol extinction is modelled using an exponential decay function with altitude. The differential transmission term is determined for two optical configurations of the MULHACEN Raman lidar. The first configuration uses VR nitrogen at ∼387 nm and water vapour at ∼408 nm, while the second uses RR at 354 nm for the molecular reference and VR water vapour at ∼408 nm. A comprehensive statistical analysis of over 300 lidar inversions from the period 2010–2016 at the Granada urban station demonstrates that neglecting $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ can lead to significant uncertainties, particularly under high aerosol loading, where values of up to 18% were observed. Up to 3 km–where most atmospheric water vapour resides–uncertainties in the retrieved water vapour mixing ratio can reach 4.0% and 12% for the first and second configurations, respectively. These uncertainties further increase to 6% and 18% up to 6 km. Based on measurements at the Granada station, including aerosols in the calculations can increase the magnitude of the differential transmission by as much as 6%.

Comparisons between lidar and SP AOD show strong agreement, with a slope of 0.93 ± 0.03, an intercept of 0.04 ± 0.01, and a determination coefficient $R^2$ of 0.80. This validates the use of SP AOD as a reliable estimate of the vertical aerosol extinction. The proposed AOD-based approach using exponential decay functions yields $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ values that closely match those derived from lidar inversions. Relative deviations between simulated and lidar-derived calculations are close to zero, with relative deviations of 0.10% (0.40%) for ∼387 nm (∼354 nm) configurations. Importantly, residuals between lidar and simulated $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculations are centred near zero with low standard deviation, indicating that including the calculation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ can convert systematic uncertainties into random uncertainties. This is particularly valuable for water vapour retrievals, as random uncertainties can be statistically treated and propagated more robustly than systematic uncertainties and can be reduced by increasing the number of measurements. Neglecting $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ may therefore introduce important uncertainty in water vapour mixing ratio retrievals.

For two different case studies corresponding to the two optical configurations of the MULHACEN Raman lidar system, calibration constants (K) for water vapour retrieval were determined with and without including the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation. In Case I, neglecting the calculation resulted in relative underestimations of 3.3% in water vapour content relative to the GNSS value. When the $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ calculation was included, this underestimation was reduced to 0.7%. In Case II, the underestimation was more pronounced, reaching 6% in water vapour content. Including the calculation reduced these differences to 0.2%, which is more critical if compared with the first configuration.

The proposed methodology enables the operational calculation of Raman lidar water vapour profiles using sun photometer data, making it suitable for long–term monitoring, climate research, and continental–scale water vapour characterisation. It is especially relevant for networks like EARLINET/ACTRIS, whose initial focus has been on aerosol studies but which can integrate water vapour as an operational product. Implementing these differential transmission calculations could improve the reliability and consistency of water vapour measurements, promoting their integration into regional and global climate assessments. Moreover, this approach reduces the dependency on the considerable effort required for lidar inversions to obtain aerosol extinction profiles, simplifying data processing and enhancing automation potential. By minimizing systematic biases associated with differential atmospheric transmission, more accurate trend detection and model validation can be performed. Nevertheless, the exact equation for the approximation of $\Delta \mathcal{T}(z,{\lambda }_{MR},{\lambda }_{H_{2}O})$ must be fitted for each station because it ultimately depends on planetary boundary layer height and typical aerosol regime of each station.