Retrieval of characteristic parameters for water vapour 1 transmittance in the development of ground based Sun-Sky 2 radiometric measurements of columnar water vapour. 3

19 Sun-sky radiometers are instruments created for aerosol study, but they can measure in the 20 water vapour absorption band allowing the estimation of columnar water vapour in clear sky 21 simultaneously with aerosol characteristics, with high temporal resolution. A new 22 methodology, is presented for estimating calibration parameters (i.e. characteristic parameters 23 of the atmospheric transmittance and solar calibration constant) directly from the sun-sky


Introduction
Water vapour columnar content is an important parameter to be estimated since it is a greenhouse component affecting the Earth climate.Many techniques were developed for measuring the water vapour amount from satellite remote sensing, in the visible, infrared or microwave spectral regions, from ground based remote sensing, i.e.GPS, sunphotometers, microwave radiometers, or from radiosondings.Sun-sky radiometers are instruments designed for the aerosol study, and many of them can also measure in the water vapour absorption band, allowing estimation of the columnar water vapour in clear sky condition, simultaneously with aerosol characteristics, with high temporal resolution up to few minutes.
Despite the limits of sunphotometry technique related to clear sky daytime conditions, the high temporal sampling and the wide distribution of these instruments all over the world make the development of methodologies for retrieving columnar water from sun-sky radiometers of great interest.The most important problem in using these instruments is the estimation of the solar calibration constant and of the a and b parameters characterizing the atmospheric transmittance in the water vapour band, b W m a e T ) ( ⋅ ⋅ − = (Bruegge et al., 1992), where m is the optical airmass and W is the columnar water vapour content.Some methods for estimation of W from sun-sky radiometers have been already developed (Halthore et al 1997, Alexandrov et al, 2009, Schmid et al, 2001).They are mainly based on the combined use of a radiative transfer code to determine the a and b parameters and of the Langley plot techniques for estimation of the solar calibration constant.Within the AERONET sun-sky radiometers network (Holben et al, 1998) a methodology for estimating W from the solar irradiance measured at wavelength of 940 nm has already implemented.Their algorithm is based on a use of a radiative transfer code (Smirnov et al., 2004) for computing T as a function of W and then estimating a and b parameters from a curve-fitting procedure.The solar calibration constant is determined by a modified Langley plot calibration performed at Mauna Loa Observatory (3400 m a.s.l).The uncertainty on its retrieval was found to be 10 times greater than the other wavelengths in the visible region, varying from 3% to 5% (T.Eck personal communication).A problem connected with these methodologies is that only one pair of (a,b) parameters is used for each kind of 940 nm interference filter, neglecting the dependence of T on the vertical profile of temperature, pressure and moisture at the various sites.This method is convenient for a network consisting of several instruments, but its correctness needs more investigations.Campanelli et al., (2010) presented a new methodology for estimating a and b parameters directly from the measurements themselves, not relying on any radiative transfer calculation and therefore reducing simulation errors and potentially containing information on seasonal changes in vertical profiles of temperature, air pressure, and moisture occurring at each measurement site.To retrieve the calibration constants from the proposed methodology some seasonal independent measurements of W (such those by radiosondes, microwave radiometers or GPS receivers) taken over a large range of solar zenith angle simultaneously with the sunsky radiometer measurements are needed.In the previous paper, data of radiosondes were used for retrieving calibration constants only in summer time, but it was also considered when such independent measurements are not available.In the latter cases, the Surface Humidity Method (SHM) was developed allowing the application of the procedure using W estimated by only measurements of surface temperature, pressure and relative humidity.
In the present paper we will improve and elaborate several points left opened in the previous paper: the study of a,b variation as a function of columnar water vapor amount by applying the methodology to an entire year dataset; the estimation of a and b retrieval errors using a Monte Carlo method; the development of a preliminary check on the quality of both sun-sky radiometer and the independent water vapour datasets; the retrieval of calibration constants using, as independent dataset, the high temporal resolution water vapour measurements from GPS receivers; the validation of the SHM examining in detail accuracy, problems and utility of this methodology.Results will be compared against measurements taken by a microwave radiometer, radiosondes and GPS receivers.

Equipment
The present methodology was applied to measurements performed during 2007 at the Chiba University (140.124 E 35.622N,34 km SE from Tokyo, Fig. 1) by the Center for Environmental Remote Sensing, Chiba University, Japan.A PREDE sun-sky radiometer model POM 02, part of Skynet network (Takamura and Nakajima, 2004; http://atmos.cr.chiba-u.ac.jp/), was used.This instrument is a scanning spectral radiometer taking measurements of solar direct and diffuse irradiance every 5 minutes at several wavelengths in the visible and near infrared regions ( 340 nm, 380 nm, 400 nm, 500 nm, 870 nm, 1020 nm) appropriately chosen for aerosol study therefore clear from gas absorption.
Measurements of direct solar irradiance taken at 940 nm are used for estimating the columnar water vapour content in clear sky conditions.Ancillary co-located measurements of pressure and relative humidity, needed for the application of the Surface Humidity Method, were provided by the Japan Meteorological Agency.
Columnar water vapour estimation from two GPS receivers stations (Shoji Y., 2013), provided by the Meteorological Research Institute of Ibaraki, Japan, were considered: n. 950225 (called GPS1 from now on) located at 35.657 N,alt.8.284 m) about 19km W from Chiba University, and n. 93025 (called GPS2 from now on) located in 35.544N,alt:50.346 m) about 10 km SW from Chiba University.
Measurements taken from a microwave radiometer (MWR) and from radiosondings (RDS) were also considered.The former (co-located with the above mentioned instruments) is a Radiometrix WVR-1100 portable water vapour passive radiometer measuring microwaves radiation from the sky at 23.8GHz and 31.4GHz.These two frequencies allow simultaneous determination of integrated liquid water and integrated vapour along a selected path.In the case of water vapour and liquid water, the atmosphere is rather translucent in the vicinity of the 22.2 Ghz water vapour resonance line, and total integrated water, water vapour and phase path delay can be derived thanks to their liner dependence on the atmospheric opacity at the measuring wavelengths.The coefficients of these linear equations are determined from bilinear regression of water vapour and inferred liquid water data derived from radiosonde observations.Radiosonde measurements were extracted from the Integrated Global Radiosonde Archive (IGRA, http://www.ncdc.noaa.gov/oa/climate/igra/)that contains quality controlled  (1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100, 70, 50, and 10 hPa), surface, tropopause and significant thermodynamic and wind levels (WMO, 1986(WMO, , 1995)).
Radiosonde estimates of the Total Precipitable Water Vapour are obtained by computing the specific humidity for each level, having valid: temperature, pressure and DPD measurements and then integrating numerically the specific humidity over the vertical using a pressure weighted numerical integration scheme.

Methodology
The direct solar irradiance measurement V [ mA] taken by the sun-sky radiometer at the 940 nm wavelength in clear sky condition is related to the solar calibration constant V 0 (extra-terrestrial current mA) at the same wavelength through the following expression where (i) m is the relative optical air mass (Kasten and Young, 1989) (Bruegge et al., 1992).Once a and b have been determined, V 0 can be estimated, and W can subsequently be calculated.
Equation (1) can be also written in the form, The aerosol optical thickness τ a is estimated at wavelength λ =940 nm, according to the wellknown Ångström formula, where wavelength λ is measured in µm, α is the so-called Ångström exponent, and β is the atmospheric turbidity parameter.Parameters α and β are determined by the regression from Eq. (3) where the spectral series of τ a are retrieved by the sun-sky radiometer measurements taken at the other visible and near infrared wavelengths 400, 500, 675, 870, and 1020 nm.
In order to find the most appropriate pair of values (a, b), the following steps are followed: i) from Eq. (2b) x-values are calculated for 30 different values of b from 0.4 to 0.7 with a step of 0.01 and each time the (x, y) squared correlation coefficient is calculated; then the maximization of the (x, y) squared correlation coefficient is used to determine the best exponent b; ii) once the optimal b exponent is retrieved, the series of x-values is computed and used in Eq. (2a) where the regression line of y versus x allows the retrieval of the coefficients a and V 0 .This modified version of Langley plot (called " type-2 modified Langley") is different from the other modified Langley method described by Halthore et al. (1997) and Schmid et al. (2001) (called " type-1 modified Langley").In fact whereas the latter determines V 0 as the intercept of the straight line obtained by fitting y versus the power term m b , in the former V 0 is retrieved by plotting y versus the product x a ⋅ where x = (m .W) b .This approach largely improves the application of the Langley methods to cases where the time patterns of W is not stable.In fact the " type-1 modified Langley" assumes that y only depends on airmass, m, and that all points have the same W.When a variability of W is recorded, the neglected dependence of y on W causes a scatter of the points and introduces calibration errors and large day-to-day changes in the retrieved calibration constants.Conversely "type-2 modified Langley" gives evidence to the dependence of y on ) ( W m ⋅ and the variability of y is explained by the real variability of the product ) ( W m ⋅ , providing a better retrieval of the intercept ( lnV 0 ) also when the time pattern of precipitable water content is not stable.In Fig. 2 type-1 and type-2 Langley plot methods were used to retrieve V 0 in two cases: stable (June 13 2007) and unstable (June 12 2007) time patterns of W as measured at Chiba by the microwave radiometer simultaneously to the sun-sky radiometer.It is clear that using type-2 method, the points are less scattered especially in the case of more unstable W time pattern.In Table 1 the retrieved V 0 values are shown.The absolute difference between the V 0 values retrieved by the two methods from 12 to 13 June is only 1.8 % if type-2 is used, whereas increases up to 4.1 % when type-1 is adopted, highlighting the better capability of type-2 in estimating V 0 during both stable and unstable W time periods.
Once parameters V 0 , a and b have been determined, the values of precipitable water content W P , can be calculated according to the equation: With respect to the previous version published in Campanelli et al. 2010, the procedure was improved in two main aspects: the use of a Monte Carlo method for the evaluation of errors affecting the a, and b retrievals and the study of their variation as a function of columnar water vapor amount by applying the methodology to an entire year dataset.Concerning the first aspect the improvement consists in: 1) A preliminary check on the quality of both sun-sky radiometer and the independent water vapour datasets (as described in Sect.4) performed before the application of the methodology .

2)
After the optimal values of a and b are found, the residual standard deviation RES σ is computed around the optimal regression line.
3) A Monte Carlo approach is used to simulate 80 fictitious bivariate samples of the pair of variables W m x ⋅ = 1 and y, each fictitious sample sharing with the true sample : the number N of data available (ii) the lower and upper bounds of x 1 (iii) the noise around the ideal straight line.
More precisely, the x 1 -data are generated by sorting N random values uniformly distributed between

MIN
x 1 and MAX x 1 , while the y-data are generated by the formula: , where the a and b (and then ln(V 0 ) ) values are the optimal values retrieved above for the given real sample, while the noise is a Gaussian noise with standard deviation coincident with RES σ .Then for each of the 80 fictitious samples a search of the optimal a and b values is carried out using the same procedure followed to find the actual optimal values (i.e., by maximizing the determination coefficient R 2 of the regression line) .In this way, a list of 80 pairs ( b a, ) are retrieved.For each of these parameters it is then possible to evaluate both the mean and the standard deviation.The coincidence of the two means a and b with their respective ideal values is a test for the goodness of the optimization procedure.This coincidence has been successfully verified in all our Monte Carlo simulations.Given that, the standard deviation ( that is the uncertainty associated to the above mean values) appears to be the best estimate of standard error to associate to each of the actual optimal values opt opt b a , , and therefore the best estimate of the uncertainty associated to the entire procedure.This evaluation is an improvement respect to the estimation obtained using a simple propagation error formula.

4)
Optimal V 0 is calculated by the linear fit of Eq.2a using the pair For what concerns the second improvement, that is the study of a,b variation as a function of columnar water vapor amount , it is evident that since b a, are supposed to depend on vertical profiles of temperature, air pressure, and moisture their "seasonal" estimation is incorrect, since seasons are only a rough subdivision of the year, marked by changes in weather, measurement environment, and hours of daylight.Therefore b a, were provided for several water vapor classes and their number their thresholds will be described in Sect 4.

Preliminary check of dataset
Simultaneous measurements of sun-sky radiometer V and independent dataset W were selected for the application of Eqs. ( 2a) and (2b).All the estimations of W within 15 minutes before and after measurements of signal V were taken, and all the values of τ A and τ R within the same intervals were selected and averaged over 30-minute time-intervals.The present method was applied in the range of solar elevation angle yielding m < 8.
A preliminary check on the quality of each dataset was performed as follows: 1. Data corresponding to τ a (940 nm) > 0.4 are rejected.
In the present study the cloud screening is performed by selecting only measurements whose RMS deviation between measured and reconstructed diffuse sky irradiance, in the wavelengths devoted to aerosol study, is lower than 8%.This criterion assured the rejection of cloud-contaminated direct and diffuse irradiance measurements, but it could not exclude the contamination of high and thin cirrus clouds.Being the maximum average value of τ a (500 nm) about 0.6 , and considering the corresponding values of Angstrom exponent, it is likely that data having τ a (940 nm) > 0.4 are contaminated by clouds, and for this reason they must be rejected, even if some good data will be probably lost .

2.
Data taken before 13:00 local time from October to May were rejected.
During these months the behaviour of y vs x appears very often not linear, as shown in Fig. 3.
In these cases two separate behaviours can be recognized generally one in the morning and one in the afternoon.This is likely related to the fact that in these months and in this time of the day (conversely to summer season) more time is needed to break the stable conditions characterizing the low atmosphere after the nocturnal cooling period.As a consequence, the vertical distribution of water vapour is anomalous respect to the profiles generally used in the in the development and/or initialization of retrieval methods (e.g.microwave radiometer, GPS, SHM) and an error can be introduced in the estimation of W. For these months we decided, as first approximation, to select only measurements initiating from 13:00 local time in order to reduce the problem to a linear behaviour.

3.
A statistical selection was applied to discard outliers with deviation greater than 2 σ from the regression line.

Parameters estimation
Because a and b are supposed to depend on the vertical distribution of the columnar water vapour and then on its total amount: i) the entire yearly independent W dataset was divided in four classes: [0-10] mm; [10][11][12][13][14][15][16][17][18][19][20] mm;  mm; [ > 40] mm; an overlap between classes of ±1 mm has been considered for the thresholds of each class; ii) the procedure was applied for each class with the aim of providing water vapour dependent a and b.The choice of a larger interval for the third class is strictly related to the need of having a great number of dataset, comparable or greater than the other three classes .
Two different independent W datasets were used for the retrieval of calibration parameters: i) The first choice is headed by the consideration that GPS is actually able to provide the more high quality estimation of W, even if a small dependence on vertical profile of temperature and water vapour needs to be corrected by an empiric relation generally retrieved from the local climatology (Shoji Y., 2013;Ortiz et al., 2011;Bevis et al., 1992).However it is not yet very common finding GPS estimations close to measurement sites, neither W from radiosondes taken over a large range of solar zenith angle, as in our case for Tateno station, where only one radiosonde launch is performed during daytime.In this case the Surface Humidity Method (SHM) can be used.
MWR in Chiba and RDS in Tateno were used for validating the results.

i) W from GPS as independent dataset
As already stated in Section 3, two stations equipped of GPS receivers are available close to Chiba University.A preliminary comparison between their results (W GPS1 , W GPS2 ) showed a difference always below 1% for all the four classes with the exception of the third class where it was found to be 2%.We decided to use W GPS2 as independent dataset for the application of the methodology, and W GPS1 for estimating the error affecting the retrieval of water vapor from sun-sky radiometer (W P ).

ii) W from SHM as independent dataset
The SHM consists in estimating W dataset using surface-level observations of moisture parameters that are much more common than those performed with radiosondes or microwave radiometers.According to Hay (1970) there is a linear dependence between precipitable water content (W SHM ) and water vapour partial pressure e 0 [hPa] at the surface, expressed by Eq. 5 where the quantity e 0 is calculated as the product of surface relative humidity f 0 by the saturation water vapour pressure E(T 0 ) [hPa], calculated as a function of surface temperature T 0 [K] according to the following Lowtran code formula (Kneizys, et al, 1983): from different daily or monthly data-sets and from varying numbers of measurements and sites (for example Hay 1970, Tuller 1977, Choudhury 1996, Liu 1986).W SHM , as defined in Eq. 5, was estimated using c 1 and c 2 coefficients taken from Yamamoto et al. (1971).They retrieved an empirical formula for the relation between W SHM and e 0 (Eq.6) using aerological measurements taken between 1950-1970 at several Japanese stations, during clear sky conditions The SHM is able to provide reliable estimation of precipitable water content when vertical humidity generally decreases as a function of height in a nearly exponential profile, but this assumption is not always verified.Undoubtedly an error in W SHM estimation can affect the validity of Eqs.( 2), but the precipitable water content amount by sun-photometric observations (W P ) can be derived accurately through Eq. ( 4), unless a and b coefficients are too far from reality, as it will be discussed in Sect.As expected the uncertainty on the determination of a and b parameters is greater for the case of SHM due probably to the lower accuracy of W SHM estimation.However in all the classes for both use of W GPS2 and W SHM the uncertainty is below 3% and 5% for b and below 9% and 14 % for a.
Looking at the water vapour dependence of the a,b and V 0 parameters in Fig. 4, is particularly noticeable that their behaviours are somehow connected since the increase of one parameter is balanced by the decrease of another.This is due to the fact that in the applied methodology of maximization these variables are not calculated independently one from the other.It implies that the slight dependence of V 0 on the water vapour class is a fictitious tendency, and therefore, at the present stage, the retrieved V 0 should be considered as an effective calibration constant whose temporal variation could not be related to a real instrumental drift.
Nevertheless its total uncertainty (estimated as the standard deviation of the values divided to their mean) resulted to be about 6% and 7% respectively when W GPS2 and W SHM are used, that is slightly largest than the maximum uncertainty retrieved by AERONET at Mauna Loa Observatory (5%).
A comparison between the two methodologies showed a general good agreement of a and b values that are always within the estimated error, with the exception of the first class where the SHM provides too low value of b and consequently an high value of a.
The behaviour of b and a as function of W is nearly parabolic with an opposite curvature.It is worthwhile recalling that the parameter a is the absorption coefficient of the water vapor band within the range 930-950 nm, weighted by both spectral curves of interference filter transmission and sensor responsivity, and that b is dependent on the intensity of the band within the spectral interval covered by sun-sky radiometer filter centered at 940 nm.The mutual correlation between W, its vertical distribution and the temperature vertical profile can affect both the parameter a ( because of the broadening of the absorption line) as well as b.
We observe that the lowest and highest W classes have a similar behavior.Such boundary classes, conversely to the other atmospheric situations, are characterized by a trapping of W due to winter inversion ( in the first one) and by the occurrence of convection (in the forth one), that favors the development of a vertical structure having one well mixed layer at the bottom and a rapid decrease upward.
In order to test such hypothesis using the available radiosonde vertical profiles, we introduced two indices to describe the W vertical distribution and having different sensitivity to the shape of the distribution.One index (P50) is the Pressure P at which is found 50% of total W. The second index is the pressure P weighted for the mixing ratio value q, (PQ) as in Eq. 7: where N is the number of vertical available measurement, taken below 100 hPa with the threshold that there are at least 16 vertical measurement to obtain a good quality radiosonding.
PQ index shows a greater sensitivity to the presence of well mixed layers respect to the P50 index, being able to discriminate (the total amount of W being equal) if water vapor is distributed within one layer or homogeneously along the entire vertical.In the former case PQ assumes values lower than the latter case.Consequently the analysis of the difference (P50 -PQ) will assume higher values when the W vertical structure will be characterize by one well mixed layer at the bottom and a rapid decrease upward.In

Water vapour estimation
Once the optimal parameters (a,b) and V 0 are estimated for each of the selected water vapour classes, a calibration table proper of the site and of the instrument under study, is made available.W P can be instantaneously calculated as in Eq.4 using this table, as soon as V (940 nm) and τ a (940 nm) measurements are performed.To retrieve the water vapour content, an iterative procedure has been set up as follows: i) for each V (940 nm) and τ a (940 nm) measurement, W P is calculated using the four set of parameters; ii) each of the four W P values falls in one class of water vapour: when at least three of them converge within the same class, the pertinent parameters to be used for the current measurement are identified.
W P was calculated using both the independent datasets from GPS (W P/GPS2 ) and SHM (W P/SHM ) and the errors affecting the retrievals (∆W P %) were estimated by a comparison against W GPS1 .The calculated absolute median percentage difference (shown in Table 2) varies from 1% to 5% for W P/GPS2 and from 1% to 11% for W P/SHM. .
The comparison between W P/GPS2 and W P/SHM (Table 3), showed a very high total correlation (0.99), and a median percentage difference varying from -0.4% (for the fourth class) up to -9% (for the third class), although always within the error ∆W P/SHM .A general underestimation by W P/SHM is observed.The most unexpected result is the small difference between the two W P estimations in the first class, where conversely the retrieved a,b parameters are very different.This topic will be deeply discussed in Section 6.
Before validating W P retrievals against MWR or RDS, we checked the goodness of these former water vapour evaluations respect to GPS (specifically W GPS2 , being closest to Chiba University where the MWR is located) that, as already stated, is actually the methodology providing the higher quality estimation of W. .We decided to correct MWR estimation by shifting W MWR values according to this formula.After the correction the disagreement between W RDS and W MWR was found to vary from 1% to 19% whereas the disagreement between W GPS2 and W MWR was found within 1% to 6%, being the higher value for the first W class.
The validation of the proposed methodology was performed by comparing W P/SHM against the corrected W MWR , W RDS and W GPS2 , whereas W P/GPS2 was compared only against the formers two.Simultaneous measurements within ± 15 minutes and ± 1 hour respectively were selected.It must be taken into account that only W RDS measurements taken at 9:00 local time can be compared with W P estimations.Scatter plots of W P/GPS2 and W P/SHM versus W MWR W RDS and W GPS2 are shown in Fig. 7c)-g) and the corresponding correlation coefficients and median percentage differences are indicated in Table 3.
W P/GPS2 and W P/SHM were found to be very well correlated with both W MWR, W RDS and W GPS2 (total correlation varying from 0.97 to 0.99).The median difference between W P/GPS2 and W MWR showed a very good agreement always within the percentage error ∆W P .The same results are found with the comparison against W RDS , with the exception of the first class were a difference of -7% was found, with a slight underestimation of W P/GPS2 respect to W RDS (0.60 mm).
The median difference between W P/SHM and W P/GPS2 showed a very good agreement always within ∆W P/SHM .The comparison against W MWR and W RDS highlighted an underestimation by W P/SHM in the second W class (-1.34 mm) and in first W class (-0.50 mm) respectively.It is worthwhile noticing that for the third class the largest disagreement was found ( -9%) showing an underestimation from W P/SHM of about 3 mm, but this class is also characterized by the greatest ∆W P (11%).
In order to validate and verify the improvements brought by the principal assumption of the proposed methodology, that is the dependence of a, b from water vapour amount, a simulation of the transmittance was performed using a radiative transfer code written by A.
Uchiyama.The code calculates the atmospheric transmittance using a correlated k-distribution method with band with 10 nm, that is a good approximation for our study.The data base of correlated-k distribution was calculated based on HITRAN data base using line-by-line code.
The code takes into account the curvature of the earth and the refraction of solar path, and  2).V 0 was calculated using the Type-2 modified Langley applied to five days having a smoothed water vapour diurnal time pattern and daily average values covering the range between 5 to 35 mm.Their mean value and standard deviation was performed to calculate V 0 s = 2.33 •10 -4 (mA) with an uncertainty of 3.5%.a s , b s and the mean V 0s value resulted to be comparable with values provided by the SHM methodology in the class having the highest water vapour content (Fig. 6).a s , b s and V 0s were therefore used in Eq. 4 to estimate the columnar water vapour W s.
The improvement taken by the hypothesis of a, b pairs dependent on W respect to the commonly used assumption of fixed a, b values, was evaluated by comparing both W S and W P/GPS2 (that is the best estimation obtained from the proposed methodology) against water vapor measured by GPS, being the most accurate retrieval of W. For this comparison W GPS1 was clearly chosen.Results (Figure 8) show that in all the W classes the agreement with W GPS1 improves when the hypothesis of variable a, b is assumed.This important outcome validates the goodness of the proposed methodology and hightlights the capability of the presented methodology of monitoring the time change of a and b values, during years, on each site and then monitoring the instrumental condition.

Discussion
Looking at the water vapour dependence of the a,b and V 0 parameters in Fig. 4 it is noticeable that the b value for class [0-10] (and therefore also a, and V 0 ) from SHM are too different (in particular too low) respect to the value retrieved using W GPS2 as independent dataset.
Nevertheless W P/SHM and W P/GPS2 for this class are in good agreement with a median difference of 3%.To explain this effect, a study of the Jacobian elements from the derivative of Eq.4 for the coefficients a and b has been performed.The analysis showed that the Jacobian for the a coefficient is approximately 3 times the one for the b coefficient, being both negative.
Therefore any sets of a and b coefficients can introduce the same error in W P determination, if the difference between b values is up to -3 times the difference between a values.When this rule is respected, two pairs of a and b can provide exactly the same W P .In our case the difference between optimal b from SHM and b from GPS dataset is about -1.5 times the difference between the corresponding a values, and this explains the good accordance between W P/SHM and W P/GPS2 .This analysis takes to the conclusion that there is a non-unique solution in the application of the SHM, unless we identify which vertical profiles of water vapour are able to provide such low b values during winter time.The problem is likely linked to the non linearity of y vs x during this season, and needs to be investigated in a next future through simulations studies.
The application of SHM needs the determination of the coefficients explaining the linearity between precipitable water content and water vapour partial pressure.In the case under study we used an empirical formula for the relation between W and e 0 obtained from a climatologic study typical of Japan, but this kind of study could not be always available for every site.One solution to this problem could be determining the proper coefficients in Eq. 6 by using already existing historical datasets of W and e 0 measured in proximity of the site under study.If this is neither possible, estimation of coefficients c 1 and c 2 can be found in the literature, from different daily or monthly data-sets and from varying numbers of measurements and sites.
With the aim of checking the error introduced if not appropriate coefficients are used for the estimation of W SHM , the Choudhury (1996) formulation was considered.Choudhury examined a data set consisting of monthly mean values of W and e 0 taken at 45 stations distributed over the entire planet, obtaining the average global values c 1 = 1.70 and c 2 = −0.1.The stations were far from water surfaces, with negligible influences due to evaporation and transport of humid air from marine regions, that are conditions not respected at all in the site under examination.Ws HM calculated by Choudhury formulation, was used in Eqs 2a and 2b for estimating the best a, b and V 0 in each water vapour class, and water vapour from the Sun-sky radiometer (W PC , where the subscript c stands for discrimination from the W P retrieved using Yamamoto formulation) was calculated using Eq.4.Linear fitting of the scatter plot between W P and W PC ( Fig. 9 a) shows an intercept value of -1.22 and a slope value of 0.92.This result indicates that even though the application of SHM can affect the validity of Eqs. 2 when completely not appropriate parameters are used for estimating W SHM , this inaccuracy introduces an error mostly consisting of a bias, positive in our case.This is also confirmed by the scatter plot of the (normalized) time derivatives of the W P and W PC time series ( ) in Fig. 9 b).In fact the optimum agreement between the two series shows that W P and W PC have the same temporal behaviour.Therefore in the case when the absolute calibration ( in terms of a, b, V 0 ) is not correct, information from the relative values of W PC and its time derivatives can be extremely valuable, being the temporal resolution of measurements high ( generally between 5 to 10 minutes).However it is strongly suggested to not use formulas of linearity between W and e 0 obtained for sites with characteristics completely different from the place under study.
Before having columnar water vapour estimations, a new user installing a sun-sky radiometer for the first time must wait to built a statistically significant dataset to the calibration table proper of the site and the instrument under study showing the water vapour dependence of the optimal parameters a,b and V 0 ,.This time is needed to collect simultaneous measurement of direct solar irradiance and pressure, temperature and relative humidity ( in the case when SHM method is used) or other independent measurements provided that they cover the entire range of variability of columnar water vapour typical of the site under study.This gap could be filled at the beginning of operations by using the method based on the simulation of transmittance, and data can be later reprocessed once the calibration table is available.
An innovative application of the presented procedure could be the possibility of providing an estimation of water vapour scale height using W P and the water vapour obtained at the ground from pressure temperature and relative humidity measurements, provided for example from the installation of sensors on the head of the PREDE sun-sky radiometer.In fact where W 0 is the water vapour density at the Earth's surface , z 0 is the scale height (km) and C is a constant taking into account the unit of measurements conversion.By inverting Eq. 8 it is possible to determine z 0 and therefore providing a sort of vertical profile from an instrument that typical retrieves only columnar properties.

Conclusions
A new methodology for determining columnar water vapour from Sun-sky radiometers measurements of direct solar irradiance at 940 nm has been introduced, based on the hypothesis that the calibration parameters characterizing the atmospheric transmittance at this wavelength, are dependent on vertical profiles of temperature, air pressure, and moisture occurring at each measurement site.To obtain calibration parameters from the proposed methodology, some seasonal independent measurements of water vapour taken over a large range of solar zenith angle simultaneously with the sun-sky radiometer measurements, are needed.In the present paper we used two independent W datasets: one estimated from GPS receivers and the other from Surface Humidity Method, a cheap procedure, easy to implement that is able to retrieve columnar W using measurements of surface temperature, pressure and relative humidity.Several aspects were developed respect to the previous paper Campanelli et al., (2010): the dependence of calibration parameter (a, b) on columnar water vapor amount for an entire year dataset; the estimation of a and b retrieval errors using a Monte Carlo method; the goodness and weakness of the SHM examining in detail accuracy, problems and utility of this methodology.
The behaviour of a and b parameters as function of W was found to be nearly parabolic with an opposite curvature.The lowest and highest W classes have similar behaviour probably because they are characterized by a W vertical structure having a well mixed layer at the bottom and a rapid decrease upward.This hypothesis was confirmed by the analysis of the available radiosonde measurements.W P obtained using GPS independent measurements, W P/GPS2 , was characterized by an error (∆W P/GPS2 ) varying from 1% up to 5% whereas W P from SHM, W P/SHM , showed an error (∆W P/SHM ) from 1% up to 11%, depending on the W classes.
The yearly time patter of W P retrieved using both the two independent W datasets, was compared against simultaneous measurements taken by a microwave radiometer, MWR, radiosonde, RDS, and GPS receivers, showing a total correlation varying from of 0.97 up to 0.99.
The accordance between W P/GPS2 and both MWR and RDS was found always within the error ∆W P/GPS2 , with the exception of the first class for RDS were a slight underestimation by W P/GPS2 (0.6 mm) was found.W P/SHM showed a good agreement with GPS retrievals, always within the uncertainty ∆W P/SHM .Its comparison respect MWR and RDS highlighted an underestimation by W P/SHM in the second W class ( -1.34 mm) and in the first W class ( -0.59 mm), respectively.
The improvement in the W P estimation brought by the assumption of a, b dependent on W was validated calculating water vapour (W S ) by using the most common procedure adopted for example by AERONET network, that consists in retrieving a and b parameters from a fitting procedure of simulated transmittance versus the product mW.W P/GPS2 and W S were compared against GPS retrievals and results showed a clear improvement using the dataset obtained by the present methodology.
Although the problems connected to the application of the Surface Humidity Method (independency of the a,b and V 0 retrievals, determination of the coefficients explaining the linearity between W and e 0 ) W P/SHM was found in good agreement with the product from different instruments.In the case when the absolute calibration ( in terms of a, b, V 0 ) resulted to be not correct, information from W P relative values and time derivatives can be anyway extremely valuable.
In conclusion we retain that the simultaneous use of the simulation method and the proposed methodology can be one solution to make the water vapour product from the sun-sky radiometer healthy, because the latter method can monitor the instrumental condition through estimation of the time change of a and b values on each site.Moreover the advantage of having simultaneous measurements of aerosol characteristics and water vapour columnar content with a high temporal resolution and obtained by using only standard surface meteorological observation for calibrating the instrument, can be of great interest to the scientific community.The present procedure will be in the next future applied to the instruments that are part of SKYNET (Takamura and Nakajima, 2004; http://atmos.cr.chibau.ac.jp/) and ESR (Campanelli et al., 2012; http://www.euroskyrad.net/)networks in which web page the open source software will be released.It will also be tested on AERONET sunsky radiometers in order to compare the two methodologies.
radiosonde and pilot balloon observations at over 1500 globally distributed stations (I.Durre et al 2006).The station closer to Chiba is Tateno ( 140.13E , 36.05 N), in the prefecture of Ibaraki about 46 km N from Chiba.The information and sampling of the radiosondings contained in the IGRA archive are, in the majority of the cases, the ones originally sent to the Global Telecommunication System (GTS) of the World Meteorological Organization (WMO).The reported variables, in the IGRA dataset, are pressure [Pa], geopotential height [m], air temperature [°C], Dew Point Depression (DPD) [°C], wind direction [°] and speed [m/s].Air temperature and DPD are reported with a 0.1°C numerical discretization.Quality assurance flags are given, for each pressure, geopotential height, and temperature value, that indicates whether the corresponding value was checked by procedures based on climatological means and standard deviations.Concerning the vertical sampling in the reported profile, in accordance with WMO guidance, radiosondes should report: standard pressure levels function of the solar zenith angle; (ii) τ a and τ R are the aerosol extinction optical thicknesses and molecular Rayleigh-scattering at 940 nm, respectively; and (iii) 940 nm as a function of m and W, with a and b constants affecting V o is obtained by evaluating the standard error on the regression line intercept ( ln(V 0 ) ) and then applying a simple propagation error formula.
weights of water vapour [ g/mol].Estimation of coefficients c 1 and c 2 can be found in the literature, 6. Calibration parameters a, b and V 0 for each W class retrieved using both W GPS2 and W SHM are in Fig 4 and Table 2.In Fig. 5 plots of the type-2 modified Langley for each of the four water vapour classes, in the case of W GPS2 , are shown.
Fig. 6 the quantity (P50 -PQ), averaged over the same four W classes analyzed in this study, is shown.It is evident that in the first and fourth class the index has the same behavior, as it happens for a and b in Fig 4, validating our hypothesis.
Figure 7 a) and b) show the scatter plot of W GPS2versus W RDS and W MWR , respectively.The disagreement with radiosonding varies from 1% to 10% (being the higher value for the first W class) with a general overestimation from RDS.Conversely the comparison against W MWR highlights a bias respect to W GPS2 , almost constant for all the classes, and expressed by the linear relationship 34 not include aerosol and cirrus clouds.The filter response function of the PREDE POM 02 was sampled; six original atmospheric models from McClatchey's (1972) (tropical, mid latitude summer, mid latitude winter, subarctic summer, subarctic winter, U.S. standard atmosphere 1962) and four modified profiles obtained by reducing the column water vapour of one tenth, were used to calculate the transmittance at 10 different hours in order to simulate a large range of path length.100 pairs of mW and simulated transmittance ( used to calculate by a fitting procedure the following parameters: a s =0.141, b s =0.626 (Table
. The grey point refers to the retrieval obtained by a fitting procedure of a simulated transmittance.

Figure 5 :
Figure 5: Type-2 modified Langley for each of the four water vapour classes.Open circles are points discarded from the quality check selection.W is from GPS2; b is the retrieved optimal value for each class, see table 2. .

Table 1 .
Retrieval of V 0 for a stable and unstable water vapor time pattern cases, using type-1 and type-2 modified Langley methods.

Table 2 .
Number of data points for each classes; optimal values of a, b and V 0 for each water vapor class and their estimated errors; estimated uncertainties of W P

Table 3 .
Correlation coefficients and median difference among W P/SHM , W P/GPS2 , W S , and measurements taken by GPS, microwave radiometer and radiosondes.