Compatibility of different measurement techniques of global solar radiation and application for long-term observations at Izaña Observatory

A 1-year inter-comparison of classical and modern radiation and sunshine duration (SD) instruments has been performed at Izaña Atmospheric Observatory (IZO) located in Tenerife (Canary Islands, Spain) starting on 17 July 2014. We compare daily global solar radiation (GSRH ) records measured with a Kipp & Zonen CM-21 pyranometer, taken in the framework of the Baseline Surface Radiation Network, with those measured with a multifilter rotating shadowband radiometer (MFRSR), a bimetallic pyranometer (PYR) and GSRH estimated from sunshine duration performed by a Campbell–Stokes sunshine recorder (CS) and a Kipp & Zonen sunshine duration sensor (CSD). Given that the BSRN GSRH records passed strict quality controls (based on principles of physical limits and comparison with the LibRadtran model), they have been used as reference in the inter-comparison study. We obtain an overall root mean square error (RMSE) of ∼ 0.9 MJm−2 (4 %) for PYR and MFRSR GSRH , 1.9 (7 %) and 1.2 MJm−2 (5 %) for CS and CSD GSRH , respectively. Factors such as temperature, relative humidity (RH) and the solar zenith angle (SZA) have been shown to moderately affect the GSRH observations. As an application of the methodology developed in this work, we have re-evaluated the GSRH data time series obtained at IZO with two PYRs between 1977 and 1991. Their high consistency and temporal stability have been proved by comparing with GSRH estimates obtained from SD observations. These results demonstrate that (1) the continuous-basis intercomparison of different GSRH techniques offers important diagnostics for identifying inconsistencies between GSRH data records, and (2) the GSRH measurements performed with classical and more simple instruments are consistent with more modern techniques and, thus, valid to recover GSRH data time series and complete worldwide distributed GSRH data. The inter-comparison and quality assessment of these different techniques have allowed us to obtain a complete and consistent long-term global solar radiation series (1977–2015) at Izaña.


Introduction
The Earth's radiation budget is essential for driving the general circulation of the atmosphere and oceans, and modulating the main conditions of the Earth's climate system.Several studies have examined the evidence of links between recent changes in climate and in the amount of daily global solar radiation (GSR H ) reaching the Earth's surface, observing a decrease in the GSR H at the surface between the 1960s and the 1990s, an effect known as dimming (see, e.g., Ohmura and Lang, 1989;Gilgen et al., 1998;Stanhill and Cohen, 2001;Wild, 2009;Wild et al., 2005), with a general decline between 4 and 6 % decade −1 over 30 years considering worldwide distributed stations.In contrast, an increase of the GSR H between 1 and 10.7 % decade −1 since the 1980s has been documented (e.g.Wild, 2012;Wild et al., 2005Wild et al., , 2007Wild et al., , 2008;;Gilgen et al., 2009), known as brightening.However, Published by Copernicus Publications on behalf of the European Geosciences Union.
for a better understanding of the global effects in the climate system, long-term GSR H data time series in representative regions are fundamental.
The first solar radiation instruments were designed in the first decade of the last century (Moll, 1913;Chaldecott, 1954;De Bruin et al., 1995;Stanhill, 1998), however, continuous GSR H observations began in the 1920s at selected locations.The longest known GSR H data time series have been measured in Stockholm (Sweden) since 1923 (Stanhill and Cohen, 2001), Wageningen (Netherlands) since 1928 (De Bruin et al., 1995) and Potsdam (Germany) since 1937 (Wild, 2015).However, regular and coordinated GSR H measurements were not well established until the framework of the International Geophysical Year in 1957and 1958(IGY, 1957/1958;Nicolet, 1982).The efforts made to increase the knowledge and measurement of the GSR H led to the necessity to centrally collect the GSR H measurements.In 1964, the World Meteorological Organization (WMO) established the World Radiation Data Centre (WRDC), which has been operating for over 50 years supported by the Main Geophysical Observatory of the Russian Federal Service for Hydrometeorology and Environmental Monitoring, and centrally collects archives and published radiometric data from several national meteorological and hydrometeorological services and other organisations (WMO, 1965).Furthermore, the increase in the comprehension of the GSR H influence in the climate system required a homogenisation in the accuracy of the measurements.Thus, several GSR H measurement networks were established around the world in the early 1990s.In 1992 the Baseline Surface Radiation Network (BSRN; http://www.bsrn.awi.de) was created (Heimo et al., 1993;Ohmura et al., 1998).The BSRN is characterised for meeting strict quality control and quality assurance protocols, and consists nowadays of approximately 60 stations in diverse climatic regions across the globe, whose radiation measurements are widely used in climate models and satellite calibration algorithms (Heimo et al., 1993;Ohmura et al., 1998).Other examples are the Atmospheric Radiation Program (ARM; Ackerman and Stokes, 2003), set up in 1992 or the Surface Radiation Budget Network (SURFRAD; Augustine et al., 2000) created in 1993.
Unfortunately, before the establishment of the cited networks, long-term GSR H data time series were very scarce and not representative for a wide variety of atmospheric conditions.In addition, throughout history, instruments of different types have been used for measuring GSR H or to obtain indirect estimations of GSR H , as sunshine duration (SD), with different accuracies, and even different radiometric scales -all this has caused inconsistencies and inhomogeneities in radiation data series, mainly before the 1950s (Fröhlich, 1991).On the other hand, new instruments with better accuracies have regularly replaced older ones, resulting in new uncertainties since in many instrument replacements there were no simultaneous measurements to assess the compatibility between old and new instruments and the corresponding changes caused in data series.
In order to complete and extend the GSR H data time series, ancillary measurements are often used to estimate GSR H .However, it is necessary to know the accuracy of those estimations by comparison with simultaneous GSR H measurements performed with modern instruments.The SD has been widely used by applying the wellknown Ångström-Prescott equation (Angstrom, 1924(Angstrom, , 1956;;Prescott, 1940) to estimate GSR H . Several authors, such as Almorox et al. (2005), Yorukoglu and Celik (2006) and García et al. (2014b), have used this method in different regions, obtaining similar results.The Robitzsch bimetallic pyranometer (also know as pyranograph or actnigraph, hereafter PYR), designed in the early 1920s (Robitzsch, 1926) and widely used until the late 1960s, measures GSR H from an equation involving the recorded area and the ambient temperature (Robitzsch, 1932).Stravisi (1986) performed a posteriori calibration of a PYR over a 3-year period obtaining hourly, daily and monthly correction factors.Later, Esteves and de Rosa (1989) proposed a correction method to improve the accuracy of daily averaged GSR H readings, reducing the error from 20 to ∼ 4 %.Maxwell et al. (1999) performed a comparison between GSR H estimations from a PYR and GSR H measurements with an Eppley PSP radiometer.They applied an automatised process to scan the PYR charts finding differences in daily GSR H values ranging between 2 and 10 % over the course of a year.
However, all these partial inter-comparisons were performed in several sites with different environmental conditions, and different instruments and time periods.In contrast, what we propose in this study is to know the performance of different instruments running in parallel in a test-bed site where the environmental conditions show a wide range of variation throughout the year.This allows us to obtain comprehensive and consistent assessments on the GSR H differences obtained with these instruments.
In this context, this work compares simultaneous groundbased GSR H measurements performed by different instruments with GSR H derived from SD and PYR area measurements in order to (1) document the traceability of main solar radiation techniques historically used, and thereby (2) assess their suitability for completing and recovering GSR H data time series, valid for climate studies.The Izaña Atmospheric Observatory (IZO) is an optimal station to carry out this quality assessment study, since solar radiation observations have been continuously performed since the early 1920s.At IZO, the SD observations started in 1917 with a Campbell-Stokes sunshine recorder (henceforth, CS recorder), which was recently replaced by a Kipp & Zonen sunshine duration sensor (henceforth, CSD recorder).The GSR H measurements started in 1977 with a PYR and it was replaced in 1992 by different instruments .Since 2005, IZO has been a station member of BRN managed by the Spanish National Radiometric Cen-tre (NRC-AEMET) and since 2009 IZO has been part of the BSRN (BSRN station #61, IZA; García et al., 2012García et al., , 2014c)).This work is divided into six sections.Section 2 describes the main characteristics of the IZO test site.Section 3 shows the technical description of ground-based instruments and the methodology used to derive the GSR H estimations from the different observations.Section 4 presents the results obtained in a 1-year intensive campaign ad hoc designed to compare different GSR H measurement techniques, and in Sect. 5 the GSR H data time series performed with two pyranometers between 1977 and 1991 at IZO is evaluated.A summary and the main conclusions of this work are given in Sect.6.
2 Site description IZO (http://izana.aemet.es) is located on the island of Tenerife (Canary Islands, Spain; at 28.3 • N, 16.5 • W, 2373 m a.s.l.), and is managed by the Izaña Atmospheric Research Centre (IARC), which forms part of the Meteorological State Agency of Spain (AEMET).
IZO is a high-mountain observatory above a strong subtropical temperature inversion layer, which acts as a natural barrier for local pollution, and it provides atmospheric measurements representative of free troposphere conditions of the North Atlantic region due to the quasi-permanent subsidence regime typical of the subtropical region (Cuevas et al., 2013;Gómez-Pelaez et al., 2013).It is optimal for solar radiation measurements given that it has about 3436 h of sunshine per year and a mean annual number of cloudfree days of 221 (61 %) between 1916 and 2015, as well as for calibration and validation activities due to a high atmospheric stability, stable ozone total column amounts, very low column water content and low aerosol concentrations.Due to these privileged measurement conditions IZO has developed a comprehensive atmospheric monitoring programme.In 1984 IZO became a member of the World Meteorological Organization (WMO) Background Atmospheric Pollution Monitoring Network (BAPMoN) and in 1989 it became a Global Atmosphere Watch (GAW) station.Also, it has actively contributed to international radiation networks and databases such as the Network for the  1).A CM-21 Kipp & Zonen pyranometer from the BSRN programme was used as reference.In the following sections the different instruments are described as well as the methodology used to derive the GSR H .

BSRN station
The GSR H measurements taken as reference in this study are those from the Izaña BSRN.Since 2009, IZO has been part of the BSRN (#61, IZA; García et al., 2012; 1).Pyranometers integrate radiation hemispherically over a horizontal surface covering a spectral range from 310 to 2800 nm (95 %).The measurements are sampled at 1 Hz, and 1 min mean values are recorded.This instrument has been calibrated recently at the World Radiation Center (WRC) at Davos, Switzerland, and it is regularly compared at IZO with a PMO6 absolute open cavity radiometer (reference instrument) designed at the Physikalisch-Meteorologisches Observatorium Davos (PMOD) following the ISO 9059:1990 (E) and ISO 9846:1993 obtaining less than 0.1 % difference between the given calibration coefficient obtained in WRC and IZO.
The expected uncertainty is ±2 % for instantaneous, hourly and daily totals (Ohmura et al., 1998;McArthur, 2005).Similar values were found by García et al. (2014c), they present a comparative study of GSR H measurements and simulations at BSRN Izaña station.The analysis was based on cloud-free days between March 2009 and August 2012 (386 days), including both aerosol-free and Saharan almost-pure mineral dust conditions.They observed agreement within 99 % and the bias (simulations-measurements) was −0.30 ± 0.24 MJm −2 (−1.1 ± 0.9 %) and RMSE of 0.38 MJm −2 (1.3 %).However, Stoffel (2005) estimated larger uncertainties for field GSR H measurements, 10 Wm −2 (±6 %).These uncertainties are associated with radiometer calibration and measurement system installation, operation and maintenance.The BSRN establishes strict quality controls for shortwave surface radiation and other measurements (Long and Dutton, 2002;Long and Shi, 2006).Applying the BSRN quality controls aforementioned to the Izaña GSR H measurements for solar zenith angles (SZA)< 90 • , García et al. (2014c) shows that the measurements largely satisfy the quality control recommended by the BSRN obtaining less than 0.1 % of the measurements out of the physically possi-ble and extremely rare limits (see Table 1 of García et al., 2014c).Also, other validation of the GSR H measurements is routinely made by comparison with the LibRadtran model (http://www.libradtran.org;Mayer and Kylling, 2005) reporting a high consistency between simulations and measurements of GSR H (see Table 5 in García et al., 2014c).

Bimetallic pyranometer
The PYR was designed by Max Robitzsch in 1926 (Table 1), and it was a popular instrument for measuring daily amount of solar irradiance from about 1932 to 1970 at meteorological stations around the world, since it is easy to install and operate and has low maintenance.The pyranometer is based on the properties of a bimetallic strip mounted under a protective glass dome, which bends in response to temperature changes.The different thermal expansion properties of the materials used to form the bimetallic strip create a simple bimetallic thermometer that responds proportionally to changes in solar irradiance (Vignola et al., 2012).By coating the strip with a black absorbent material and exposing it to the Sun's energy, the strip bends in proportion to the change in temperature, which is proportional to the intensity of the solar radiation.The signal is instantaneously recorded in a strip card, being the area under the radiation curve proportional to the incoming solar radiation.The measurement uncertainties introduced by a variety of undesirable response characteristics result in daily solar irradiance within ±20 % for the best-performing instruments (Coulson, 1975;Garg and Garg, 1993).For this reason this instrument was classified by WMO as a third-class pyranometer (WMO, 1965).
In order to convert the area obtained from the pyranometer strip card to radiation units the following equation is used: where K is the instrument constant, S is the value of the integrated surface (area) in cm 2 and t is the averaged temperature in • C from temperature observations at 07:13:18 UTC.Q is expressed in cal cm −2 min −1 .The constant K is dependent on each instrument; unfortunately this constant is unknown for the instrument used in this work.Therefore, we have estimated the constant K from Eq. ( 1) considering that Q is the theoretical GSR H .In particular, we have calculated a K value for each month considering only cloud-free days.These days were selected using Long and Ackerman's method (Long and Ackerman, 2000).Hereafter, the GSR H estimations with this method are denoted PYR-1 GSR H . Also, we have used a second approach, the Esteves method (Esteves and de Rosa, 1989), to determine the GSR H performed with the pyranometer (hereafter, PYR-2 GSR H ).This method defines a monthly factor (F m ) that is determined for a given set of cloud-free days of the same month as the sum of ratios between the theoretically calculated daily GSR H value and the area of the register for the same day, divided by the number of days in the set: where H ci is the theoretical GSR H for the clear day i, S i is the area inscribed under the registered curve for the day i and n m is the number of clear days of the month considered.
Once the value of the monthly factor is obtained, it is possible to calculate the daily GSR H for any day of the month as follows: where S is the area inscribed under the register's curve for a given day of the month (idem to Eq. 1).GSR H is expressed in MJm −2 .For the two methods, two variables are needed to be calculated: the theoretical GSR H (Q in PYR-1 and H ci in PYR-2) and the area S. Firstly, the theoretical GSR H was calculated using the LibRadtran model (see Sect. 3.1).To simulate the GSR H values, we used daily aerosol optical depth (AOD) data from AERONET (Holben et al., 1998) and precipitable water vapour (PWV) observations from Vaisala meteorological radiosondes (Miloshevich et al., 2009), launched at Güimar (WMO station #60018) at 11:15 UTC (Romero Campos et al., 2011).Further details about the LibRadtran simulations can be found in García et al. (2014c).
The area, S, was calculated by using a semiautomatic method ad hoc developed at IZO to digitalise the data recorded on the pyranometer strip cards.Figure 1 shows the flow chart illustrating the methodology used on 10 August 2014, as an example.The processing is performed using the MATLAB Image Processing Toolbox (Russ, 1994).Firstly, each strip card is scanned with a resolution of 300 ppi and saved to 24 bit RGB TIF (Tagged Image Format).The image is trimmed in order to only include the radiation curve in the image, discarding possible spurious notations and the signals on the strip card.In a few cases, minor editing of the trimmed image is needed to clean it of spots and defects.A median filter is applied to the trimmed image and two background areas are selected (Fig. 1.3), normally in the upper corners of the image.For each channel (R, G and B) a statistical analysis is performed on the background areas in order to establish the thresholds needed to define the pixels that belong to the curve.Once the curve pixels are identified we extract the coordinates in pixel units and the curve is smoothed to detect the sunrise and sunset points using the Savitzky-Golay algorithm (Gorry, 1990;Fig. 1.6).With the pixel coordinates of the sunrise and sunset we construct a straight line connecting both points.Finally the area, S, in pixels 2 is obtained by calculating the area bounded between the radiation curve and the base line.Taking into account that the resolution of the image is 300 ppi, dividing the area S in pixels 2 by 13 950, the area in cm 2 is obtained.We have also computed the instantaneous GSR H values from the pyranometer radiation curve.The conversion of the horizontal axis from pixels to time is straightforward by taking into account the relation between the distance in pixels and time from sunrise to sunset throughout the day.The vertical conversion requires more effort due to the clear dependence of the calibration factors on month and hour, as pointed to in previous works (Albrecht, 1954;Stravisi, 1986).Thus, assessment of calibration factor dependence on the mentioned variables is required.In order to convert from cm to Wm −2 the next relation is proposed: where Y is obtained in the image processing method already discussed and C are coefficients computed with the LibRadtran model which has demonstrated to be a very useful measurement quality control tool (García et al., 2014c).This study is performed for a sample of cloud-free days in May and June (2015), a month in which there is a maximum SZA range.

Sunshine recorders: CS and CSD
The CS recorder was invented by J. F. Campbell in 1853 and modified to its usual shape by G. G. Stokes in 1879 (Table 1).
The CS focuses the direct solar beam through a glass sphere, mounted concentrically in a section of a spherical bowl, on a burn card located under the sphere.The card is provided with a time indication, which makes it possible to determine the SD from the length of burn when the card is removed from the instrument at the end of the day (Painter, 1981;WMO, 1996;Sanchez-Lorenzo et al., 2013;García et al., 2014a, b).The errors of this recorder are mainly due to the dependence of burning initiation on card's temperature and humidity as well as to the overburning effect, especially in the case of broken clouds (Kerr and Tabony, 2004).
The CSD recorder (Table 1) was designed with the objective of automating the process of measuring SD for meteorological services and has replaced the traditional CS recorders.This instrument is formed by three detectors: the sensor at the front measures global radiation, while the other two detectors, at the middle and at the rear, are partially shaded covering 1/3 of the sky, measuring diffuse radiation (Kipp and Zonen, 2003).Furthermore, the three detectors have exactly the same spectral and angular characteristics, making the process of calibration easy.The direct solar radiation can be computed from the values of global and diffuse solar radiation, and the sunshine duration is determined according to the latest WMO definition (WMO, 1996).
In this work, the method used to determine the GSR H estimated from SD records was developed by Ångström (Angstrom, 1924(Angstrom, , 1956) ) and later modified by Prescott and Rietveld (Prescott, 1940) and applies the following equation: where H and H o are the daily GSR H and the daily extraterrestrial GSR H (MJm −2 day −1 ), respectively, on a horizontal surface, n and N d are the number of hours measured by the SD recorder and the maximum daily SD, respectively, and a and b are coefficients to be determined by using linear regression.This method has successfully been validated at IZO by García et al. (2014a, b), obtaining an expected uncertainty of 9.2 % for CS recorder and 5.5 % for CSD recorder when comparing BSRN GSR H .Following Eq. ( 5) and details given in García et al. (2014a, b), the a and b coefficients at IZO are given in Table 2 for each SD recorder.To account for the variability introduced by the presence of clouds, fog, etc., these coefficients were estimated as functions of the fraction of clear sky (FCS), which is defined as the ratio between the maximum daily sunshine duration N d and SD performed with CS recorder (n):

MFRSR radiometer
The MFRSR (Table 1) was developed in the early 1990s.This radiometer uses independent interference-filter photodiode detectors and the automated rotating shadowband technique to make spectrally resolved measurements at seven channels with six wavelength passbands of 10 nm FWHM centred near 415, 500, 610, 665, 862 and 940 nm, and an unfiltered channel used to obtain broadband solar irradiance estimates from 300 to 1100 nm.Each measurement sequence is routinely repeated every 15 s and the recorded instantaneous signals are subsequently averaged over a 1 min time interval to increase the signal-to-noise ratio (Harrison et al., 1994).
In this work, the MFRSR was re-calibrated using the Li-bRadtran model obtaining a new calibration factor calculated as the mean of ratio between the GSR H measurements and simulations in a time interval of 5 min before and after the solar noon (WMO, 1996), only considering cloud-free days (160 days, 44 %).The obtained calibration factor was 0.438 ± 0.008.
Although the MFRSR provides the GSR H directly, the spectral range covered by this instrument (300-1100 nm) is significantly smaller than that of the reference instrument used, the CM-21 pyranometer (280-2600 nm).For this reason, we have calculated a radiation correction using the Li-bRadtran model as a linear relation between the GSR H in the MFRSR and the CM-21 spectral ranges.We have modelled the same days selected to calculate the calibration factor in both spectral ranges with the same input values, obtaining a very good agreement (correlation coefficient, R, of 0.995) and the following fit equation: GSR 280−2600 = (1.239± 0.006)GSR 300−1100 + (0.4247 ± 0.014).( 7

Results
In this section we present the comparison of the daily GSR H values obtained with the different instruments and techniques using BSRN as reference and we also perform an analysis of the GSR H bias as a function of the season, solar irradiance, temperature, relative humidity (RH), FCS and AOD.
In addition, we present the comparison of the instantaneous GSR H values measured with PYR and MFRSR with respect to BSRN.

Comparison of daily GSR H data
Applying the different methodologies discussed in the previous section, we have computed the daily GSR H values from the different instruments and techniques between 17 July 2014 and 12 July 2015 at IZO.The daily GSR was calcu-lated according to the following equation: where GSR E (t) is the instantaneous value of the solar irradiance (Wm −2 ) at time (t), computed from sunrise (sr) to sunset (ss).
In general, as expected, the best agreements with BSRN GSR H (with R = 0.999) are found for those instruments that directly measure solar radiation (not estimations) and have been calibrated with the LibRadtran model, since the model has been widely validated against the BSRN measurements at IZO (García et al., 2014b).Although the poorest scores are those provided by GSR H estimations from CS measurements (R = 0.978), it is worth emphasising that the correlation coefficient is always greater than 0.97 and the slope ranges between 0.92 and 1.04 for all cases.The correla- tion between PYR-1 GSR H and PYR-2 GSR H is very good, with R = 0.997.Results of comparison of the instruments are shown in Fig. 2.
In order to quantify the uncertainty of the different techniques with respect to BSRN GSR H , we have calculated the absolute difference or bias (XXX GSR H -BSRN GSR H , in MJm −2 ) and the relative difference ((XXX GSR H -BSRN GSR H )/BSRN GSR H , in %), where XXX means PYR, CS, CSD and MFRSR.As a summary, Table 3 lists the metrics used to quantify these measurements (MB: median bias; SD: standard deviation; and RMSE: root mean square error).In general, GSR H is underestimated, except for GSR H MFRSR.The results obtained for PYR-1 and PYR-2 GSR H are very similar with MB of −0.5 MJm −2 (−2 %), while the precision, given by the RMSE, is of 0.9 MJm −2 (3.5 %) for both cases.The greatest difference is obtained for the CS GSR H , with MB of −1.14 MJm −2 (−4.5 %) and RMSE of 1.70 MJm −2 (6.7 %).The MFRSR GSR H is slightly overestimated with a MB of 0.28 MJm −2 (1.1 %) and RMSE of 0.54 MJm −2 (2.1 %).Since the results obtained by PYR 1 and PYR-2 are similar, hereafter we will discuss only the results of the PYR-1 -renamed as PYR.
As aforementioned, some of the instruments and methods analysed are sensitive to different factors and atmospheric conditions.We have analyse the GSR H differences with respect to BSRN GSR H as a function of the solar irradiance (Fig. 3a), the average temperature and RH (Fig. 3b and c, respectively), the FCS and AOD.
The PYR GSR H estimations do not show dependence on solar irradiance (Fig. 3a) with a MB almost constant around −0.5 MJ −2 .The PYR GSR H shows a slight dependence on temperature with the lowest MB in the 10-20 • C range (Fig. 3b) which agrees with the seasonal variation, considering that the highest median temperatures are found in summer and autumn at IZO.This behaviour was also observed by Stravisi (1986) and Esteves and de Rosa (1989) and is attributed to the response of the bimetallic strips to temperature changes (see Sect. 3.2), since the optimum operating temperature is around 15 • C. Figure 3c shows that there is not a clear dependence on RH.
The CS GSR H estimations show a clear dependence on solar irradiance, with higher bias for higher BSRN GSR H values.However, there is no temperature dependence, with an almost constant bias value around −1 MJ −2 (Fig. 3b), and a slight dependence on RH, with higher bias for the lower RH values (Fig. 3c).This might be explained by the fact that a card requires more solar irradiance to be burnt under wetter conditions, as frequently occurs in early morning (Bider, 1958), while it burns more easily under warm and dry conditions (Kerr and Tabony, 2004;Painter, 1981).Thus, overburning of the card in spring and summer seasons is to be expected.
The CSD GSR H estimation, unlike CS GSR H , shows much lower bias, with an opposite annual cycle to that found for CS GSR H , with higher GSR H bias in winter and autumn, and lower in spring and summer.Although a slight dependence on solar irradiance is observed, this is much lower than in the case of CS GSR H (Fig. 3a).There is a slight dependence on temperature, with higher negative bias for lower temperatures (Fig. 3b) and a MB close to zero for temperatures > 20 • C. The manufacturer declares a temperature dependence of 0.1 %/ • C ( Kipp and Zonen, 2003).There is no dependence on RH (Fig. 3c).
Finally, the MFRSR is the instrument that shows the best performance, with a bias close to zero through the whole year, and the lowest scatter (Fig. 3).We observe, in general, an overestimation in MFRSR GSR H , unlike the rest of the compared instruments.The MFRSR GSR H has a clear positive dependence on solar irradiance (Fig. 3a).There is no temperature dependence for temperatures > 15 • C, (Fig. 3c) and a slight dependence for lower temperatures.The MFRSR is thermally controlled at around 40 • C. Thus, when the difference between ambient temperature and 40 • C is very large, the heating system is continuously working at 100 % effort, creating in former MFRSR instrument versions an electro- magnetic frequency interference with the solar irradiance measurements, leading to higher measurement inaccuracies (Harrison et al., 1994;Hodges and Michalsky, 2011).The MFRSR GSR H measurements show a slight positive dependence on RH (Fig. 3c).
We have also studied the differences with respect to FCS and AOD (not shown here).No dependence on FCS was found, although it should be noted that 85 % of the days (N: 232 days) showed FCS > 75 % while only 1 % (N: 4 days) showed FCS < 25 %.Concerning the dependence on AOD only background conditions (AOD < 0.10) and dust conditions (AOD ≥ 0.10) have been considered based on García et al. (2014b).No dependence on AOD is found, although we must highlight the fact that 87 % of the days (N: 231 days) showed AOD < 0.10 and 13 % of the days (N: 33 days) showed AOD ≥ 0.10.The GSR H measurements most affected by AOD were those obtained with the CSD, showing negative bias for pristine skies and positive bias for dust conditions.

Diurnal GSR E comparison
This section presents the diurnal comparison of the instantaneous GSR E values from PYR and MFRSR (the only instru- ments that have this information) with those of the BSRN, used again as reference, applying the methodology discussed in Sect.3.2.We have considered 27 cloud-free days between May and June 2015, grouping them into intervals of SZA from 10 to 60 • (Table 4, Fig. 6).The results show that the GSR E MB decreases with the SZA for PYR GSR E , while  MB is quite stable for MFRSR GSR E values (Table 4 and Fig. 6).In both cases the SD increases notably with SZA.

1977-1991 Time series
As an application of the previous analysis, we have applied the methodology explained in Sect.4.1 to the measurements performed at IZO between 1977 and 1991 with two bimetallic pyranometers: (1) #290604 (1977)(1978)(1979)(1980)(1981)(1982), and (2) #250585 (1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991).The theoretical GSR H needed to obtain the instrument constant (K in Eq. 1, PYR-1 GSR H ) in these periods was calculated using the LibRadtran model, where the main model inputs were the AOD and PWV.Given that there are no available AOD measurements, we have used AOD derived from artificial neural networks (ANNs) applying the methodology developed by García et al. (2016).The PWV between 1977 and 1979 was taken from the NCEP/NCAR reanalysis (Kalnay et al., 1996) and since 1980 we have used the PWV obtained from meteorological radiosondes (Miloshevich et al., 2009), and specifically those launched from Santa Cruz de Tenerife station (WMO#60020).Once the instrument constants were evaluated we calculated the GSR H from the PYR measurements.By combining the re-evaluated GSR H data time series between 1977 and 1991 from PYR measurements and the GSR H series measured with different pyranometers from 1992 to 2015, obtained in García et al. (2014b), we finally completed the 38-year GSR H data time series at IZO, depicted in Fig. 5.
In order to analyse the suitability of the PYR measurements to fill gaps in long-term GSR H data time series, we have studied the temporal stability of the PYR data time series at IZO between 1977 and 1991 by comparing two independent GSR H data series.During that period, the only available GSR H data at IZO are those estimations from CS records, presented in García et al. (2014b).That work presented the re-evaluation of the GSR H data time series between 1933 and 1991 from SD measurements performed with a CS by using the Ångström-Prescott method, and documented also its high quality and temporal consistency.Hereafter, we refer to that data time series as CS-AP GSR H to distinguish it from the CS GSR H obtained in the current study.
In general, there is a good agreement between PYR GSR H and CS-AP GSR H (Fig. 6a), being both GSR H data time series very consistent with each other, with correlation coefficient of 0.86.To check the temporal stability of PYR GSR H observations, we examine possible drifts and discontinuities in the MB data time series (PYR GSR H -CS-AP GSR H ), considering a drift as the linear trend of the annual MB, while the change-points (change in the median of the bias time series) are analysed by using a robust rank order change-point test (Lanzante, 1996).We found that there are no significant drifts in the MB time series and no change-points are found at 99 % level of confidence.To complete this analysis, we ex-amined the temporal evolution of the annual transfer function between PYR GSR H and CS-AP GSR H , shown in Fig. 6b.This is calculated as the annual slope and intercept considering all the PYR GSR H and CS-AP GSR H observations for each year.As for MB, neither temporal drifts nor significant discontinuities were detected at 99 % confidence level for any of the least-squares-fit parameters.
The slope values are higher than 0.88 for all the years, with the median value of 0.92, coincident with that obtained between PYR GSR H and CS GSR H for the 2014-2015 period (Fig. 2).The intercept in the period 1977-1991 is higher than in the 2014-2015 period, while the R values are of about 0.95 in the whole period, and 0.98 in the 2014-2015 period.This improvement is likely due to the cleaning and fitting of the instrument before being restarted for the 2014-2015 intercomparison.

Summary and conclusions
A 1-year-long comparison (17 July 2014to 12 July 2015) between GSR H measurements performed with old and modern radiation and sunshine duration instruments has been performed at IZO.The daily GSR H values measured with a bimetallic pyranometer (PYR) and multifilter rotating shadowband radiometer (MFRSR), and GSR H estimated from SD performed by a Campbell-Stokes sunshine recorder (CS) and a sunshine duration sensor (CSD) have been compared with respect to GSR H from a BSRN CM21 pyranometer.We have also compared instantaneous values of GSR H performed with PYR and MFRSR for different SZA with respect to GSR H BSRN.
Assuming BSRN GSR H as reference, the measured or estimated GSR H values show median biases of 2 and 1 % for PYR and MFRSR GSR H , respectively, and of 5 and 2 % for CS and CSD GSR H , respectively.These results, as expected, show that the instruments that directly measure GSR H , such as the PYR and MFRSR, show lower MB and lower scatter than the ones that estimate the GSR H , such as the CS and CSD recorders.Moreover, MB values for each instrument are within their corresponding uncertainty, agreeing with results obtained by other authors (Coulson, 1975;García et al., 2014b;McArthur, 2005).The comparison of the daily GSR H values from PYR and MFRSR showed a good agreement with BSRN GSR H , obtaining a RMSE of ∼ 0.9 (3 %) and ∼ 0.5 MJm −2 (2 %) for PYR and MFRSR GSR H , respectively, and ∼ 1.7 (7 %) and ∼ 1.1 MJm −2 (4 %) for CS and CSD GSR H , respectively.It is worth highlighting the fact that the biases for PYR found in this study are lower than those reported by others authors.For example, Coulson (1975) and Garg and Garg (1993) obtained uncertainties between 10 and 20 % reduced up to 4-5 % by Esteves and de Rosa (1989) and Soulayman and Daudé (1995).These results, obtained with simultaneous observations under the same environmental conditions, provide information about expected GSR H uncertainties from historical instruments useful for assessing long-term GSR H data series constructed from both classical and modern instruments.
The GSR E from PYR and MFRSR were compared with GSR E BSRN on a daily basis at different solar zenith angles.The RMSE increases from 3 to 7 %, and from 0 to 3 % for PYR and MFRSR GSR E , respectively, when the SZA increases.
The methodology developed in this work has allowed us to obtain a unique re-evaluated and quality-assured GSR H data time series from 1977 to 1991 at IZO using GSR H data from two bimetallic pyranometers.The consistency in the obtained daily GSR H PYR is good with respect to the daily GSR H estimates from SD measurements with a correlation coefficient of 0.86.
These results demonstrate that (1) the continuous intercomparison of different techniques used to measure GSR H constitutes an important diagnostics tool for identifying inconsistencies between GSR H data records, and (2) GSR H measurements performed with classical and simple instruments are fundamental for filling gaps in long-term GSR H data series after accurate data screening, calibration and correction by using ad hoc reanalysis techniques based on transference radiative models.

Data availability
The daily GSRH data time series for the period 1977-2015 at IZO are available at http://izana.aemet.es/dataseries/GSR.
Competing interests.The authors declare that they have no conflict of interest.
Acknowledgements.This study forms part of the activities carried out within WMO CIMO Testbed for Aerosols and Water Vapour Remote Sensing Instruments at Izaña Observatory.The authors are grateful to the IZO team and especially all observers who have worked in the past 70 years taking bimetallic pyranometer measurements at IZO as well as all members of security staff who have helped make this work possible.
Edited by: O. Dubovik Reviewed by: three anonymous referees

Figure 1 .
Figure 1.Flow chart of the method used to extract the areas from the strip cards of the bimetallic pyranometer.

Figure 3 .
Figure 3. Box plot of bias (PYR: black; CS: red; CSD: green; and MFRSR: blue) versus (a) corresponding BSRN GSR H measurements (MJm −2 ), (b) temperature in • C and (c) relative humidity (RH) in % between 17 July 2014 and 12 July 2015 at IZO. Lower and upper boundaries for each box are the 25th and 75th percentiles, the solid line is the median value, and hyphens are the maximum and minimum values.N indicates the number the measurements in each season.

Figure 4 .
Figure 4. Box plot of the bias (PYR: black; MFRSR: blue) with respect to BSRN GSR E in Wm −2 at different ranges of SZA between May and June 2015 (27 cloud-free days) at IZO. Box plots are defined as in Fig. 3.

Figure 5 .
Figure 5. Daily GSR H data time series between 1977 and 2015 at IZO.The green and blue dots correspond to the measurements performed with PYR #290609 and #250585, respectively, between 1977 and 1991, and the grey dots represent the measurements performed with different pyranometers (CM-5, CM-11 and CM-21) between 1992 and 2015.

Figure 6 .
Figure 6.(a) Time series of the median bias between PYR GSR H and CS-AP GSR H in MJm −2 , and (b) time series of the slope (black line), correlation coefficient (red line; R) and intercept (blue line) for each year.The black, red and blue dashed lines indicate the slope, R and intercept obtained in the period 2014-2015, respectively (see Fig. 3).The green dashed lines represent the instrument change in 1983.

Table 2 .
García et al., 2014b)nts obtained between 1992 and 2000 for the CS recorder and between 2001 and 2013 for the CSD recorder, as a function of FCS (fraction of clear sky, %) and % of days used (seeGarcía et al., 2014b).SEM is standard error of the mean.

Table 3 .
Statistics of the comparison between daily GSR H from PYR (Method 1 and 2), CS recorder, CSD recorder and MFRSR with daily BSRN GSR H (XXX GSR H -BSRN GSR H ) at IZO between 17 July 2014 and 12 July 2015 (in MJm −2 ).MB, median bias; SD, standard deviation; MAB, mean absolute bias; and RMSE, root mean square error; statistics percentages in parentheses.

Table 4 .
Statistics for the bias between instantaneous GSR E performed with PYR and MFRSR with respect to BSRN GSR E (XXX GSR E -BSRN GSR E ) at IZO between May and June 2015 (27 cloud-free days) (in Wm −2 ) for different range of SZA.MB, median bias; SD, standard deviation; and RMSE, root mean square error.Statistics for the relative bias are in parentheses.