Water Vapor Retrievals from Spectral Direct Irradiance Measured with an EKO MS-711 Spectroradiometer—Intercomparison with Other Techniques

: Precipitable water vapor retrievals are of major importance for assessing and understanding atmospheric radiative balance and solar radiation resources. On that basis, this study presents the ﬁrst PWV values measured with a novel EKO MS-711 grating spectroradiometer from direct normal irradiance in the spectral range between 930 and 960 nm at the Izaña Observatory (IZO, Spain) between April and December 2019. The expanded uncertainty of PWV (U PWV ) was theoretically evaluated using the Monte-Carlo method, obtaining an averaged value of 0.37 ± 0.11 mm. The estimated uncertainty presents a clear dependence on PWV. For PWV ≤ 5 mm (62% of the data), the mean U PWV is 0.31 ± 0.07 mm, while for PWV > 5 mm (38% of the data) is 0.47 ± 0.08 mm. In addition, the EKO PWV retrievals were comprehensively compared against the PWV measurements from several reference techniques available at IZO, including meteorological radiosondes, Global Navigation Satellite System (GNSS), CIMEL-AERONET sun photometer and Fourier Transform Infrared spectrometry (FTIR). The EKO PWV values closely align with the above mentioned different techniques, providing a mean bias and standard deviation of − 0.30 ± 0.89 mm, 0.02 ± 0.68 mm, − 0.57 ± 0.68 mm, and 0.33 ± 0.59 mm, with respect to the RS92, GNSS, FTIR and CIMEL-AERONET, respectively. According to the theoretical analysis, MB decreases when comparing values for PWV > 5 mm, leading to a PWV MB between − 0.45 mm (EKO vs. FTIR), and 0.11 mm (EKO vs. CIMEL-AERONET). These results conﬁrm that the EKO MS-711 spectroradiometer is precise enough to provide reliable PWV data on a routine basis and, as a result, can complement existing ground-based PWV observations. The implementation of PWV measurements in a spectroradiometer increases the capabilities of these types of instruments to simultaneously obtain key parameters used in certain applications such as monitoring solar power plants performance.


Introduction
In recent decades, greenhouse gases (GHG) have become of major interest as main drivers of climate change [1]. The most important GHG is water vapor, since it efficiently interacts with both solar and terrestrial radiation, and its concentration in the atmosphere is driven by temperature, with a strong positive feedback in the Earth's climate system [1]. Hence, it has been defined as an "essential climate variable" by the World Meteorological Organization (WMO) Global Climate Observing System (GCOS) program [2], with its observation being critical for characterizing the Earth's climate system and its changes. Spatial and temporal changes in the distribution of the water vapor content have been pointed out as an indicator of climate change [3]. An increase in the tropospheric content of water vapor has been reported [4].
Apart from its key role in climate change, water vapor is involved in many other processes in the atmosphere; e.g., the Earth's energy budget, the hydrological cycle, and in cloud and aerosol formation [5]. In addition, water vapor is one of the greatest sources of error in space-based remote sensing and communications, leading to errors in satellite pointing and attenuation on their signals [6]. In this context, precise measurements atmospheric water vapor content measurements with high temporal frequency and spatial resolution are mandatory.
To quantify the atmospheric water vapor content, the integrated precipitable water vapor (PWV) is widely used, which is defined as the total amount of water vapor present in a vertical atmospheric column from the observation site up to the top of atmosphere. The techniques to measure the PWV are very diverse and include airborne, space, and groundbased instruments using different approaches (passive or active sensors, different spectral regions or observation geometries). Focusing on ground-based measurement techniques, meteorological radiosondes provide precise humidity vertical profiles, and then obtaining integrated PWV values that have the disadvantage of corresponding to a relatively long time interval (60-90 min). Moreover, the number of PWV observations obtained each day with this technique is relatively low, since although there are 1300 stations spread all over the world (https://www.wmo.int/pages/prog/www/OSY/Gos-components.html# upper; last access 1 December 2020), only one or two launches are made at these stations per day (most of them only one), and the cost of each sounding is relatively high. Very accurate PWV values are retrieved from microwave radiometer profilers (MWP), which measure the radiation emitted by the atmosphere in the microwave spectral range. They have the advantage of measuring under any weather condition [7], but the inconvenient of its high cost and maintenance requirements. A less precise approach, but allowing PWV data to be provided under any condition as well, is the use of the Global Navigation Satellite System (GNSS) satellites to retrieve the water vapor content by means of the delay suffered by the signal received at the ground-receptors. By means of a post-processing procedure, it is possible to convert this delay in PWV, allowing a very high spatial coverage and temporal frequency [8].
Solar radiation measurements have been successfully used for the monitoring of the Earth's gaseous composition since last century. Among existing solar measurement techniques, sun photometry is a well-developed tool for water vapor observations (e.g., [9][10][11][12][13][14][15][16][17][18] and references therein). Presently, the filter radiometers (Sun and Lunar photometers) initially devised to measure aerosols, may provide PWV retrievals by measuring the direct sunlight in spectral channels centered in the water vapor absorption bands, usually at 719, 817 and 946 nm [19]. The benefits of sun photometers (lower cost compared with other instruments and easy installation) lead them to be chosen as reference instruments to establish several worldwide networks: AErosol RObotic NETwork (AERONET; [20,21]), Global Atmospheric Watching Precision Filter Radiometer network (GAW-PFR; [22]) and Skynet Radiometers network (SKYNET; [23,24]).
More sophisticated PWV solar devices are spectrometers operating in the infrared domain such as those based on the Fourier Transform Infrared (FTIR) technique, which provide precise water vapor profiles and integrated PWV amounts by analyzing the measured high-resolution solar absorption spectra (e.g., [25][26][27]). In the visible and near infrared spectral range, PWV data can be also retrieved from direct irradiance measured by spectroradiometers, but the works on this matter are scarce (e.g., [28][29][30]). The first studies date from the beginning of the 20th century when a method to derive the PWV by spectroscopy was described [31] and later applied at Mount Wilson, California to retrieve the PWV [32]. Since then, several retrieval methodologies have been proposed to retrieve PWV from solar absorption spectra, including monochromatic approaches, spectral windows, differential optical absorption spectroscopy (DOAS), and iterative nonlinear fitting models (e.g., [10][11][12][13]29,30]).
Given the diversity of measurement techniques and approaches available in the literature, numerous comparison studies have been carried out in the last decades (e.g., [12][13][14]18,25,[33][34][35]). These studies show differences ranging between 0.40% and −25.4%. A summary of the PWV intercomparison results obtained in this study is shown in Section 6.
The main advantage of a spectroradiometer is that, with a single instrument, a set of parameters can be accurately measured (spectral irradiance, integrated irradiance, aerosol optical depth, irradiance attenuation by clouds, and PWV) which are key input parameters of certain applications, such as determining the performance of solar plants. In this context, this paper presents for the first time the PWV retrievals of a novel grating spectroradiometer, the EKO MS-711. The PWV data were determined from the spectral direct normal irradiance (DNI) measurements in the spectral range between 930 and 960 nm taken at the Izaña Observatory (IZO). The quality assessment of the new EKO PWV retrievals was carried out using the Monte-Carlo method and was complemented with the comparison to reference techniques available at IZO between 1 April and 31 December 2019: radiosonde, GNSS, CIMEL-AERONET sun photometer and FTIR. The work is divided in 7 sections: Section 3 describes the main characteristics of the Izaña station, and the instruments and techniques used in this work. Section 4 and 5 details the approach used to determine PWV and its associated uncertainty estimated from Monte-Carlo method, respectively, while Section 6 shows the comparison to the reference PWV observations. Finally, a summary and the main conclusions are given in Section 7.

Site Description
The observations used in this work were performed at the Izaña Observatory (IZO) that is part of the Izaña Atmospheric Research Center (IARC) from the State Meteorological Agency of Spain (AEMET). This sub-tropical observatory is located on Tenerife island (The Canary Islands, Spain; 28.3 • N, 16 [36]. More details of the station and the measurement programs can be found in [37]. Given the altitude and location of the IZO, it is normally above a temperature inversion layer, avoiding possible mixture processes with local pollution from low levels of the island. Consequently, if offers excellent conditions to perform in situ and remote sensing atmospheric measurements of aerosols and trace gases under free troposphere conditions. With respect PWV, IZO is characterized by very dry conditions. Ref. [35] studied the TCCON FTIR and CIMEL-AERONET PWV time series between 2007 and 2019 at IZO, observing that the wettest months occur in summer (maximum PWV in August) and the driest months correspond to winter (minimum PWV in February), with a total average PWV of 4.8 ± 2.8 mm from TCCON FTIR, and 3.8 ± 1.4 mm from CIMEL-AERONET.

EKO MS-711 Spectroradiometer
Within the framework of the Izaña Observatory WMO-CIMO testbed activities, an EKO MS-711 grating spectroradiometer was deployed and set in a DNI measurement mode. The MS-711 spectroradiometer model was originally devised to measure global solar spectral radiation. By attaching a collimation tube to the sensor head, the instrument Field of View (FOV) can be narrowed to a 5 • , and, when assembled on a solar tracker, it can measure the spectral DNI. The instrument measures in the 300-1100 nm spectral range with a bin width of 0.4 nm and a bandpass of nominally <7 nm (defined as the full width at half maximum (FWHM)) [38,39]. The instrument was installed on an EKO sun-tracker STR-21G-S2. With this configuration, the solar DNI spectrum is sampled with a temporal resolution of 1 min, with the integration time of each measurement varying between 10 ms to 5 s, depending on the light intensity. For more details refer to [40]. The PWV content in the atmosphere column is determined from the measured vertical profiles of temperature, relative humidity and atmospheric pressure acquired with the meteorological radiosondes [6]. The obtained PWV values have been post-corrected for radiation dependence (daytime radiosondes) and temperature following et al. [41]. The precision is 5% for total column water vapor in the lower and middle troposphere, increasing to between 10% and 20% in the upper troposphere for very dry conditions [6,25]. • Global Navigation Satellite System (GNSS) Using GNSS receivers, it is possible to determine the PWV content due to the delay of the signal caused by water vapor molecules [8,42]. The principle of the PWV retrieval lies in the decomposition of the indicated signal delay received from the GNSS receivers in components such as the ionosphere, the clock, the geometric and the neutral atmosphere's delays. Ref. [8] developed the technique to evaluate the atmospheric delay caused by refractivity effects, which occurs mainly in the troposphere, and split it into two components, the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD). The ZHD is caused by dry air plus the non-dipole contribution of water vapor and can be accurately modelled from the surface pressure, altitude and latitude [8]. The PWV is calculated from the ZWD by the knowledge of a weighted mean temperature of the atmosphere, as defined in [8].

Ancillary Instrumentation
In this work, PWV retrievals from a Leica GRX 1200GG pro GPS/GLONASS receiver, installed at IZO since 2008, have been used. This receiver belongs to the Spanish National Geographic Institute (IGN), and forms part of the Network of European Meteorological Services (EUMETNET) GPS water vapor program (EGVAP) [43]. The GNSS PWV retrievals came from precise orbits and are available every hour (GNSS-P; [41]). Their expected precision is of 20% for PWV < 3.5 mm and about 10% above this value [25].

• CIMEL Sun Photometer
Since 2003, IZO has formed part of AERONET that is a federation of ground-based remote sensing aerosol networks [20]. The reference instrument of AERONET is the CIMEL Electronique CE318-T multiband sun photometer [44] that performs measurements of spectral sun and moon irradiance and sky radiances, at 340, 380, 440, 500, 675, 870 and 1020 nm with a FOV of 1.3 • [45]. An additional channel centered at 940 nm (with a FWHM of 10 nm) is used to retrieve PWV, following the procedure described in [21], with an uncertainty of 10%. Ref. [25] found that the uncertainty for this instrument depends on humidity conditions, with 7% for humid conditions and 25% for very dry conditions (PWV ≤ 2 mm). In this work, we use AERONET Version 3.0 Level 1.5 data (http://aeronet.gsfc.nasa.gov, last access: 10 December 2020). • FTIR Within the IZO's atmospheric research activities, the FTIR solar measurements started in 1999 in the framework of a collaboration between the AEMET and the KIT (Karlsruhe Institute of Technology) [46], contributing to NDACC and TCCON networks since 1999 and 2007, respectively. In this work, the FTIR PWV data retrieved within the TCCON activities are used [47]. These data are derived with the 2014 version of the GGG processing software [26] by analyzing direct solar absorption spectra measured between 4500 and 6500 cm −1 (corresponding to wavelengths λ between 1500 and 2200 nm) at a spectral resolution of 0.02 cm −1 . The latter represents a resolution power λ/∆λ at 5000 cm −1 of about 2.5 × 10 5 [25]. The aperture diameter is limited to 0.5 mm, leading to a narrow spectrometer's FOV angle of only 0.07º. Several scans are co-added to increase the signal-to-noise ratio and, thus, the sampling frequency of TCCON FTIR PWV data is about 2 min at IZO. Using collocated meteorological radiosoundings at globally distributed TCCON sites, a scale factor of 1.0180 ± 0.0040 has been determined for the TCCON FTIR water vapor columnar products (XH 2 O) [26].
Recently, at IZO, Ref. [35] further studies corroborated these results, obtaining a mean bias of −0.006 mm (−1.3%) when comparing TCCON FTIR PWV values to those obtained from meteorological radiosondes launched on Tenerife island and processed according to the GRUAN scheme (http://www.gruan.org, last access: 28 November 2020). The Almansa's study also found a relative difference of +3.4% respect to collocated CIMEL-AERONET data (FTIR-CIMEL). Further details of the FTIR program at IZO are given [25,37,46].

PWV Retrieval Method
The basis of the methodology used in this work lies in the Beer-Lambert-Bouguer law [33,48]. In spectral regions of strong variation of the spectral molecular absorption, as the water vapor absorption band, this law can be modified as follows: where DNI(λ) is the direct normal irradiance at wavelength λ measured by the spectroradiometer, DN I o (λ) is the top of atmosphere irradiance corrected for the Sun-Earth distance, τ(λ) is the optical depth, m is the optical air mass and T w (λ) is the water vapor transmittance. Considering that the most important components for τ(λ) in the considered absorption band are the Rayleigh optical depth and the aerosols (Figure 1), this equation can be written as follows: where τ R (λ) is the molecular Rayleigh scattering [49] that depends on the station pressure, m R is the Rayleigh air mass calculated to according to [50], τ a (λ) is the extinction aerosol optical depth (AOD) and m a is the aerosol air mass calculated according to [51]. The dependence of T w (λ) with PWV was studied by [10,52], obtaining the following dependence: where m w is the water vapor air mass that is approximately equal to m a , PWV is the columnar water vapor content, and a and b are constants. Ref. [12] found that these constants depend on wavelength position, the shape and width of the filter functions, the atmospheric pressure and temperature at the station as well as the vertical distribution of water vapor. Substituting Equations (2) and (3) into Equation (1), we obtain the following expression: Therefore, the PWV can be determined from the following equation: The PWV has been determined in the spectral range between 930 and 960 nm (Figure 1), considering the integrated band, with a modification of Equation (3), as follows: where ∆λ = λ 2 − λ 1 : λ 1 and λ 2 are the wavelength limits (930 and 960 nm, respectively).
In preliminary analysis of this study, the approach of [30] was used, obtaining PWV retrievals very similar (not shown here) to those obtained with the approximation indicated below. Since the moderate resolution of the EKO MS-711 spectroradiometer (10 nm) in the spectral range in which the PWV is retrieved, we decided to choose a more robust method in which the whole spectral range in which the water vapor absorbs significantly is used, to minimize the signal-to-noise ratio. The PWV can be derived from as follows: The coefficients a and b have been determined from Equation (3) by a fitting plot in logarithm scale of the ln(ln(1/T w (∆λ))) against ln(m w PWV) resulting in a line with the slope equals to b, and the intercept equals to ln(a) [10] (Figure 2). To do so, the T w (∆λ) has been simulated using the MODTRAN 6 model [53] multiple runs varying the solar zenith angle (SZA) between 0º and 90º with steps of 1º and water vapor from 0 to 40 mm with steps of 2 mm steps at the site altitude of 2373 m a.s.l., and considering the standard midlatitude summer atmosphere model [54]. The modelled spectra were convolved with the slit function of the EKO MS-711. DN I o (λ) was determined using the modified Langley-plot technique [14,40] for each wavelength between 930 and 960 nm (Figure 2). The τ a (λ) has been determined from the DNI measurements performed with the EKO MS-711 following the methodology described in [40]. In that work the AOD values were evaluated at the 340, 380, 440, 500, 675 and 870 nm wavelengths, thus the spectral AOD between 930 and 960 nm has been derived by applying the Ångström formula [55]. It is important to note that the FOV of the EKO spectroradiometer is 5º, whereby the AOD retrievals have been post-corrected of circumsolar radiation following the procedure presented in [40]. The effect of the large EKO MS-711 FOV on the measured solar absorption spectra was studied in [40] by the circumsolar ratio (CR) defined as: where CSR is the circumsolar radiation coming from within the instrument FOV, and DN I SUN is the direct normal irradiance coming from the Sun disc. The CR definition allow us to provide the relation between the DN I SUN and the DN I measured with the EKO: To study the possible effect of the EKO FOV on the PWV retrievals, DN I SUN and CSR have been simulated varying the PWV between 0 and 20 mm for AOD of 0.05, 0.10 and 0.20 in the spectral range between 930 and 960 nm. The simulations were done by using the LibRadtran model [56,57], following the procedure explained in [40]. The results are shown in Figure 3 for a SZA of 30º. As can be seen in the figure, CR is under 0.8% for all the simulations done, independently of the wavelength, AOD or PWV values. The Figure 3 also points out that the most influent parameter on the CR is the AOD, as in [40]. However, for the AOD values typically found at IZO (0.05 and even lower) [37], the CR impact drops to less than 0.2%. Finally, the PWV content has no appreciable effect on the CR and, therefore, no corrections due to FOV impact on the DNI measurements have been applied on the EKO PWV retrievals.

PWV Uncertainty Analysis
This section presents the uncertainty estimation of the EKO PWV retrievals due to the propagation of uncertainties associated with the parameters and measurements involved in the PWV retrieval presented in the previous section. To do so, the standard and internationally recognized methods proposed by the BIPM (Bureau International des Poids et Mesures), as the so-called GUM guides (Guide to the Expression of Uncertainty in Measurement, [58][59][60][61]) have been followed to evaluate the uncertainty. In particular, we have adopted the widely used Monte-Carlo method (MCM), as explained in [59], and only a short outline with the most important concepts used in this procedure is given in this paper. Refer to the cited guides for a detailed description.
The associated uncertainty of a measure or variable X is stated by a distribution function, determining the probability that its value is less than or equal to any value x. The derivative of this distribution function is the probability density function (PDF), which provides the probability element, i.e., the probability that x < X < x + dx [58]. In this work we are assuming that all variables considered in the analysis are uncorrelated and for all the measured ones the uncertainty is of type-B ([58], 4.3). This implies the available information about their uncertainties is only a bound p such that the probability of the value X to be in the range [X − p, X + p] is one, and zero out of this interval ([59], 4.3.7), defining therefore a rectangular PDF distribution.
Considering a model f (e.g., an equation) such that for input parameters X an output variable Y can be obtained, in the MCM, M random input values are assigned to the inputs X, obeying the distribution defined by their PDF, and therefore obtaining M output estimators of Y (y r ) through the model f. Thus, in summary, the MCM propagates the input PDFs to obtain a discrete representation of the output PDF. The latter can be characterized by the standard uncertainty (u, [58], 7), which is the uncertainty of a measurement expressed as a standard deviation, estimated by the following expression: The number of runs (M) depends on the level of confidence desired for the uncertainty of Y. In the GUM framework, this coverage probability represents the probability that the value of Y is within a range of values with a stated probability. In the MCM, it is equivalent to the expanded uncertainty in the propagation of uncertainties obtained by the way defined in [58]. In the current work M has been set equal to 10 6 to reach a 95% of coverage probability or level of confidence ([59], 7.2). The Python module MetroloPy (https://pypi.org/project/metrolopy/; last access: 16 October 2020) has been used to perform the MCM calculations.
To estimate the standard uncertainty of PWV amounts, the MCM has been applied considering that our model f is given by Equation (7), where the input parameters (X) considered are listed in Table 1. This table also displays the uncertainty values used, their references, and the PDF assigned to each of them. Please note that previously the MCM has been also used to evaluate the uncertainty associated with the coefficients a and b by using Equation (3), and the DN I o (λ) by considering the Modified Langley-Plot technique, since no information about their uncertainty were available. The PDFs for a, b and DN I o (λ) have been assumed as normal distributions, once evaluated with the MCM, while the rest of the variables as rectangular distributions since they are obtained from references. The MCM PDFs obtained for the input variables a, b, and DN I o (λ) fits very well to normal distributions, which was corroborated by applying the Anderson-Darling test [64]. In this test, the data sample is tested against a specific distribution (normal, lognormal, etc.), being its null hypothesis that the sample follows the distribution. This test was also applied also to the MCM results of each PWV retrieval, obtaining a high agreement of the PDF to a normal distribution.
This fact allows us to associate an expanded uncertainty (U PWV ) to each PWV retrieval as two times the standard uncertainty (u) obtained from the MCM application. It is important to ensure these preliminary assumptions since if the PDF of the MCM result for each retrieval were not close to the normal distribution, large differences could be introduced in the estimated MCM coverage interval (JCGM [59], 5.11.6.d and 9.2.4).
Once the PDFs for the input parameters are stated Equation (7) is run M times obtaining M estimators of each PWV. After sorting in ascending order ( [60], 7) and applying Equation (10) to these M estimators, the standard uncertainty and the coverage interval are obtained. The Figure 4 shows a schema of the PVW standard uncertainty determination. The temporal evolution of the EKO PWV amounts between April and December 2019 and the relative U PWV is shown in Figure 5. The relative U PWV ranges between about 3% (0.13 mm) and 34% (0.81 mm) with a mean and standard deviation of 9.8 ± 4.1% (0.37 ± 0.11 mm), respectively. From this figure the dependence of the U PWV on the PWV content ( Figure 6) becomes evident: higher PWV amounts provide lower relative U PWV values. For PWV ≤ 5 mm (62% of the data), the mean of the U PWV is 11.9 ± 3.9% (0.31 ± 0.07 mm), while for PWV> 5 mm (38% of the data) is 6.5 ± 1.1% (0.47 ± 0.08 mm). Also, U PWV shows a strong dependence with AOD ( Figure 6), again the higher AOD values lead to lower relative U PWV . For AOD ≤ 0.10 (95% of the data), the mean of U PWV is 10.1 ± 4.1% and for AOD > 0.10 (5% of the data) the mean decreases until 5.9 ± 1.4%.  Please note that the highest uncertainties in PWV are obtained under pristine conditions (AOD < 0.03). For hazy conditions (AOD > 0.1), always caused by the presence of Saharan dust at the Izaña Observatory, the uncertainty in PWV is within 10%. The explanation for these results is found in the fact that the observations of PWV and mineral dust are not independent. The Saharan dust-laden air masses (higher values of AOD) are also associated with higher PWV compared to pristine free troposphere conditions [65][66][67] in which the observed PWV (≤ 2 mm) is close to the instrumental detection limit. Therefore, the percentage uncertainty of PWV under moderate-high AOD is much lower than that obtained in extremely dry conditions. The MCM analysis allows the relative contribution of each uncertainty source to be determined and, therefore, the potential improvements to be identified. The most important contributors to the uncertainty are the DNI and AOD, with a mean contribution of 49% and 42%, respectively. The DNI o (λ) contributes a mean of 9.2% and the input of the rest of parameters is less than 0.01%. Therefore, more precise direct solar measurements and AOD estimations, will yield more precise EKO retrievals.

Results: Comparison to PWV Observations
In this section, we present the comparison between the EKO PWV retrievals determined in the spectral range between 930 and 960 nm and other techniques available (RS92, GNSS, CIMEL-AERONET and FTIR) at Izaña Observatory between April and December 2019.
The comparisons were done considering the different temporal resolution of each technique available at IZO. To avoid redundant data, each EKO measurement is only paired once to the reference observations, considering the closest measurement in time for each temporal window. For RS92, as the reference time, we choose the time at a half of the observation (a sonde takes approximately one hour between its launch and its burst in the upper troposphere-lower stratosphere). Therefore, we have taken the EKO PWV values 30 min after the launch time, which is the expected time for the radiosonde to reach the IZO altitude (2.4 km) [25]. The GNSS PWV values have a temporal frequency of 1 h, thereby the closest GNSS record within 5 min around each EKO measurement have been paired. Finally, for FTIR and CIMEL-AERONET the temporal window was reduced to 2 min around each EKO measurement. This approach produced quasi-synchronous and non-duplicated comparison sets of 187 data for RS92/EKO, 1987 data for GNSS/EKO, 3765 data for FTIR/EKO and 21,909 data for CIMEL-AERONET/EKO. Figure 7 summarizes the comparison of the paired PWV observations.
To quantify the difference of the EKO PWV respect to the PWV of different techniques, we have calculated the PWV bias (EKO -XXX, in mm) and relative PWV differences ((EKO -XXX)/XXX), in %). As a summary, Table 2 lists the metrics considered to quantify the differences, for which the following expressions were used: where XXX represents the PWV of the reference techniques RS92, GNSS, FTIR and CIMEL-AERONET.
The comparisons show a good agreement between EKO PWV and the different techniques with Pearson correlation coefficients (R) higher than 0.950 (Table 2; Figure 7). The 86% of the differences between EKO and GNSS and CIMEL-AERONET PWVs are within the range ±1 mm (Figure 7), decreasing to 75% and 70%, for RS92 and FTIR, respectively.  The smallest differences respect to EKO PWV are obtained for the CIMEL-AERONET PWV with a RMSE, STD and R of 0.67 mm (12.8%), 0.59 mm (13.6%) and 0.983, respectively, while the greatest differences are found respect to RS92 PWV with a RMSE, STD and R of 0.94 mm (18.5%), 0.89 mm (25.3%) and 0.956, respectively. However, it should be noted that the best slope obtained (Figure 7) corresponds to the fit between RS92 versus EKO PWV with a value of 0.983. The difference is low between EKO PWV and CIMEL-AERONET PWV because both instruments measure directly the sunlight and the technique applied in both instruments is similar to determine the PWV, while that the RS92 is an in situ technique and therefore air mass differences may result [25].
It is important to note that the number of paired EKO/CIMEL-AERONET PWV measurements (N 21909) is much higher than for the rest of the techniques. On the other hand, when analyzing the bias among techniques, the least differences are found respect to GNSS PWV with MB of 0.02 (3.3%), while the largest differences are respect to CIMEL-AERONET and FTIR PWV with MB of 0.33 mm (9.5%) and −0.57 mm (14.5%), respectively. These systematic differences can be attributed the possible calibration inaccuracy of the AERONET PWV product [35], while the differences found respect to FTIR PWV are likely due to spectroscopy inconsistences between the different spectral ranges in which the PWV is determined (930-960 nm for EKO PWV and 1500-2200 nm for FTIR PWV). The bias found respect to different techniques are within the uncertainty determined with the Monte-Carlo method (Section 5) with values ranging between 0.13 and 0.81 mm.
The PWV relative differences are found to be dependent on the PWV content, such that they considerably decrease when comparing EKO PWV values higher than 5 mm (Table 2).
These results are within the same order of magnitude as those found by other authors at IZO (Table 3), i.e., ref. [25] found differences of −25.4 ± 12.7% between FTIR and CIMEL-AERONET PWV, −5.36 ± 19.5% between FTIR and GPS PWV, and −3.33 ± 15.5% between FTIR and RS92 PWV in the period 2005 and 2009. Ref. [34] compared the PWV obtained with Lunar CIMEL photometer at IZO between July and August 2011, and reported differences of 1.8 and 2.5 mm respect to GNSS and RS92 PWV, respectively. Ref. [35] found relative differences of 9.1% between FTIR and ZEN-R52 radiometer PWV and 17.1% between CIMEL-AERONET and ZEN PWV from August 2017 to June 2018. The results obtained in other PWV retrieval comparisons are also shown in Table 3, with similar differences than those obtained with the EKO MS-711.

Summary and Conclusions
In this study, we determined the first PWV retrievals performed with an EKO MS-711 grating spectroradiometer from normal direct irradiance in the spectral range between 930 and 960 nm at IZO, in the period April-December 2019 (Number of spectra: 30826). The EKO PWV was compared with PWV retrievals routinely obtained by four different techniques at IZO. The PWV values are in the 0.08-24.1 mm range, with most values (90% of data) below 8.6 mm.
The uncertainty estimation (U PWV ) of the EKO PWV retrievals was determined using the Monte-Carlo method, considering the propagation of uncertainties associated with the parameters and measurements involved in the PWV retrieval. The mean of UPWV is 0.37 ± 0.11 mm. The uncertainty presents a clear dependence with the PWV content. For PWV ≤ 5 mm (62% of the data), the mean is 0.31 mm, while for PWV > 5 mm (38% of the data) is 0.47 mm. These results are within the estimated range of precision established by different authors with other techniques.
The most important contributors to the PWV uncertainty are the DNI and AOD, with a mean contribution of 49% and 42%, respectively. The DN I o (λ) contributes a mean of 9.2% and m w with 0.013%. The rest of the parameters contribute less than 0.001%.
When comparing the EKO PWV with the rest of techniques available at IZO, our retrievals show a consistent agreement with them, with R > 0.950 for all comparisons. The smallest differences were obtained against the CIMEL-AERONET PWV with a RMSE and STD of 0.67 mm and 0.59 mm, respectively, while the most significant differences were found against the RS92 PWV with a RMSE and STD of 0.94 mm and 0.89 mm, respectively. Concerning the bias, we find the least differences between GNSS and EKO PWV with MB of 0.02, and the largest differences against CIMEL-AERONET and FTIR techniques with MB of 0.33 mm and −0.57 mm, respectively. These differences decrease considerably for atmosphere conditions with PWV ≥ 5 mm, obtaining MB values ranging from −0.45 mm (against FTIR PWV) to 0.11 mm (against CIMEL-AERONET PWV). These results confirm that the methodology applied in this work to obtain the PWV is consistent with the different techniques available at IZO, with at least the same order of error. This work, together with [40], demonstrate the capability of the EKO MS-711 spectroradiometer to provide high-quality AOD and PWV measurements, with excellent results when compared against well characterized instruments with established records for reliability. Considering that the EKO MS-711 spectroradiometer was designed for global spectral radiation measurements, it is remarkable that with a technical adaptation, accurate and satisfactory PWV results were possible on its DN I(λ) measurement configuration. The main advantage of a spectroradiometer is that with a single instrument, a set of parameters can be accurately measured (spectral irradiance, integrated irradiance, aerosol optical depth, irradiance attenuation by clouds, and PWV) which are key input parameters of certain applications, such as determining the performance of solar power plants.

Institutional Review Board Statement:
Not applicable for studies not involving human or animals.

Informed Consent Statement:
Not applicable for studies not involving human.

Data Availability Statement:
The EKO MS-711 data might be available upon request from EKO Instruments and Izaña WMO-CIMO test bed. The CIMEL-AERONET data from the Izaña Atmospheric Observatory ("Izana") are available on the AERONET. The TCCON Izaña FTIR PWV data are available at TCCON Caltech DATA archive at https://data.caltech.edu. The GNSS and RS92 data used are available from the authors upon request.