Reconstruction of global solar radiation time series from 1933 to 2013 at the Izaña Atmospheric Observatory

This paper presents the reconstruction of the 80year time series of daily global solar radiation (GSR) at the subtropical high-mountain Izaña Atmospheric Observatory (IZO) located in Tenerife (The Canary Islands, Spain). For this purpose, we combine GSR estimates from sunshine duration (SD) data using the Ångström–Prescott method over the 1933/1991 period, and GSR observations directly performed by pyranometers between 1992 and 2013. Since GSR measurements have been used as a reference, a strict quality control has been applied based on principles of physical limits and comparison with LibRadtran model. By comparing with high quality GSR measurements, the precision and consistency over time of GSR estimations from SD data have been successfully documented. We obtain an overall root mean square error (RMSE) of 9.2 % and an agreement between the variances of GSR estimations and GSR measurements within 92 %. Nonetheless, this agreement significantly increases when the GSR estimation is done considering different daily fractions of clear sky (FCS). In that case, RMSE is reduced by half, to about 4.5 %, when considering percentages of FCS > 40 % (∼ 90 % of days in the testing period). Furthermore, we prove that the GSR estimations can monitor the GSR anomalies in consistency with GSR measurements and, then, can be suitable for reconstructing solar radiation time series. The reconstructed IZO GSR time series between 1933 and 2013 confirms change points and periods of increases/decreases of solar radiation at Earth’s surface observed at a global scale, such as the early brightening, dimming and brightening. This fact supports the consistency of the IZO GSR time series presented in this work, which may be a reference for solar radiation studies in the subtropical North Atlantic region.


Introduction
Solar radiation controls the energy radiative balance in the Earth and, thus, our weather and climate.For this reason, its study has been one of the main objectives of the research community during the last decades.Recently, the focus is on evaluating the long-term trends of solar radiation reaching the Earth's surface (global solar radiation, GSR) as well as on identifying the variability driven by climate change (Stanhill and Cohen, 2001;Sanroma et al., 2010;Wild, 2009).Observational evidences of changes on GSR trends have already been reported at a global scale.A decrease of the solar radiation at the surface has been observed between the 1960s and the 1990s, an effect known as dimming, with a general decline between 4 and 6 % over 30 years considering worldwide distributed stations (e.g.Ohmura and Lang, 1989;Gilgen et al., 1998;Stanhill and Cohen, 2001;Liepert, 2002;Pinker et al., 2005;Wild et al., 2005;Sanchez-Lorenzo et al., 2007;Wild, 2009).In contrast, since the 1980s a partial recovery has been documented, with an increase of the solar radiation known as brightening (Wild et al., 2005(Wild et al., , 2007(Wild et al., , 2008(Wild et al., , 2012;;Gilgen et al., 2009).Trends between +1.0 and +10.7 % decade −1 have been reported from the 1980s onwards (Wild, 2009, and references herein).
The causes of the dimming/brightening phenomena are not well understood yet (IPCC, 2007;Wild et al., 2012), Published by Copernicus Publications on behalf of the European Geosciences Union.
but some authors associate them with changes in global cloud cover, atmospheric aerosols content (natural or anthropogenic), and to the role of direct and indirect aerosol effects (Stanhill and Cohen, 2001;Ramanathan et al., 2001;Wild et al., 2005;Wild, 2009).The relative importance of these factors may differ depending on region and pollution level (Wild, 2009).In this context, some authors point out that the dimming may only be a local effect, associated with urban environments (e.g.Alpert and Kishcha, 2008), or that tendencies over land and over ocean can differ in sign and in magnitude (e.g.Pinker et al., 2005).For a better understanding of these global effects and to reduce the uncertainties that still remain, long-term GSR time series in regions representative of background signals are fundamental.
Reliable solar radiation studies need high-quality, longterm and worldwide distributed GSR measurements.Although 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), regular and coordinated GSR observations were not well established until the 1950s within the framework of the International Geophysical Year, with expected accuracies of 5 % (Nicolet, 1982a, b).Currently, different radiation programmes and international networks (WMO, 1984) ensure the high quality and consistency of GSR records by using very precise instrumentation (uncertainty of 2 %, WMO, 2008).In order to complete gaps in these GSR time series, to correct erroneous records, or to extend them over time, GSR estimations from other climate variables, such as sunshine duration (onwards, SD), cloud cover or visibility, are very valuable.For this purpose, the most extended approach is to estimate the GSR from SD data (Iziomon and Mayer, 2002;Sivamadhavi and Selvaraj, 2012), since it combines long-term records availability (SD measurements started in the 19th century, Butler and Hoskin, 1987;Pallé and Butler, 2001;Sanchez-Lorenzo and Wild, 2012), simplicity and reliable results.The relation between SD and GSR observations have been described by using several mathematical relations, including linear, the so-called Ångström-Prescott relation (Angstrom, 1924;Prescott, 1940), cubic (Samuel, 1991), logarithmic (Ampratwum and Dorvlo, 1999) and exponential (Almorox and Hontoria, 2004) relations.In all of these equations, a set of coefficients are calculated in a simultaneous period of SD and GSR measurements.The comparison of the aforementioned approaches conclude that the use of complex relations instead of the simple linear relation proposed by Ångström-Prescott does not significantly improve the GSR estimates (Almorox and Hontoria, 2004;Yorukoglu and Celik, 2006).
In this context, the goal of this work is to reconstruct the GSR time series between 1933 and 2013 at the Izaña Atmospheric Observatory (IZO) by using the Ångström-Prescott method on SD measurements (1933( /1991( ), and GSR observations performed by pyranometers (1992( /2013)).For this purpose, this work is divided into six sections.Section 2 describes the different instruments and measurements used (radiation and SD data) as well as the main characteristics of IZO.Section 3 shows the recalibration of measured GSR at IZO between 1992 and 2005, when no strict quality controls were applied on GSR measurements, by using the Li-bRadtran model (onwards, LibRadtran).Section 4 explains the method applied to estimate the GSR from SD records and documents the precision and consistency of our GSR estimates, while the reconstruction of the whole GSR time series from 1933 to 2013 is addressed in Sect. 5. Finally, a summary and conclusions are given in Sect.6.

Site description, measurements and tools
This study has been performed at the high-mountain Izaña Atmospheric Observatory (IZO, http://izana.aemet.es),located in Tenerife (The Canary Islands, Spain; 28.3 • N, 16.5 • W, 2373 m a.s.l.).This station is managed by the Izaña Atmospheric Research Center (IARC) from the Meteorological State Agency of Spain (AEMET).For almost three decades, the IARC has aimed to monitor atmospheric constituents that are capable of forcing changes in the climate of the Earth, through modification of the atmospheric radiative environment (greenhouse gases and aerosols), and those that may cause depletion of the global ozone layer.IZO has been part of the WMO-GAW programme (World Meteorological Organization-Global Atmospheric Watch) since 1984 and part of the NDACC (Network for the Detection of Atmospheric Composition Change) since 2001.Also, it actively contributes to international aerosols and radiation networks such as AERONET (Aerosol Robotic Network) and BSRN (Baseline Surface Radiation Network).IZO is a suitable site for in situ and remote sensing observations, optimal for calibration and validation activities due to a high atmospheric stability, high frequency of clean and pristine skies, a stable total column ozone, very low column water content and low aerosol concentrations.IZO provides atmospheric measurements representative of free troposphere conditions of the subtropical North Atlantic region due to the quasipermanent subsidence regime typical of the subtropical region (Cuevas et al., 2013;Gomez-Pelaez et al., 2013, and references therein).

Global solar radiation data
Between 1992 and 2013 the GSR measurements have been performed with different pyranometers at IZO (Table 1).
-Before 2005, the measurements were taken with two instruments (Kipp & Zonen CM-5 and CM-11), but they were not subject to strict quality controls.According to the manufacturer, daily uncertainties of ±5 and ±2 % can be expected, respectively (http://www.kippzonen.com).
-Since 2009 IZO has been part of the BSRN (BSRN station #61, IZA).In 2004 the BSRN was designated as the global baseline network for surface radiation for the Global Climate Observing System (GCOS) with the main objective of providing measurements of the best quality for shortwave and long-wave surface radiation fluxes, with high temporal resolution.For this purpose, the BSRN establishes exigent quality controls for its products (Long and Dutton, 2002;Long and Shi, 2006).
Applying the BSRN quality controls aforementioned to the Izaña GSR measurements for solar zenith angles (SZA) < 90 • , R. D. García et al. (2012) showed that the IZO measurements largely satisfy the quality controls recommended by the BSRN.The estimated uncertainty for the instantaneous global irradiances indicated by BSRN in 1997 is ±5 W m −2 (±2 %) (Ohmura et al., 1998;McArthur, 2005).These values account for calibration uncertainties and were estimated from standard deviation of the calibration coefficients.
Since 2009 the BSRN and NRC programmes have coexisted at IZO, keeping their own instrumentation and data evaluation procedures.They show an excellent consistency: Pearson correlation of 0.98 and RMSE of the daily differences in the common period 2009/2013 of 0.34 W m −2 (1.1 %), within the instrumental uncertainty.Note that the different pyranometers acquire GSR records on a 1 min basis.However, in this work, we use daily GSR values, calculated by integrating the 1 min measured GSR from sunrise to sunset (García et al., 2014).
To complete the gaps observed in the long-term GSR time series at IZO, we have used the GSR measurements taken at the Teide Observatory (OT, http://www.iac.es)managed by the Instituto de Astrofísica de Canarias (IAC).OT (28.3 • N, 16.5 • W) is only 1.3 km far away from IZO and at the same altitude (2371 m a.s.l.).The GSR measurements at OT were performed with a Silicon Cell Pyranometer (SCP, Model 3120), with an expected uncertainty of ±3 % for daily values (see http://www.allweatherinc.com/).The periods in which we used these measurements were April-August 2000, January 2001-July 2002and September 2003-July 2005.Nonetheless, a difference < 0.5 % between the GSR in the range 250-1500 nm and in the range 310-2800 nm is observed.

Radiation transfer model and input parameters
LibRadtran is a complete free software package containing a suite of tools for radiative transfer calculations in the Earth's atmosphere (freely available from http://www.libradtran.org)(Mayer and Kylling, 2005).
The GSR and DSR (direct solar radiation) estimations for IZO were computed by using this model as described in García et al. (2014).The radiative transfer calculation is addressed by a multi-stream discrete ordinates algorithm (Stamnes et al., 1988) (DISORT for SZA ≤ 70 • and SDIS-ORT for SZA > 70 • ) and the input model parameters (atmosphere gas composition, surface albedo, aerosol optical properties, . . . ) are directly measured at IZO.
The most significant changes with regard to García et al. (2014) are related to aerosol optical depth (AOD) and total precipitable water vapour (PWV) observations.Since 1994 the AOD values have been measured with different sun-photometers.In this work we have used AOD data from the following instruments (in brackets the expected uncertainty is given for each instrument): a PMOD/Rocket between 1994 and 1996 (±0.013 at 368 nm, Díaz et al., 2000), a MFRSR (Multifilter Rotating Shadowband Radiometer) between 1996and 2001(±0.03, Pedro Miguel Romero, personal communication, 2014), a PFR (Precision Filter Radiometer) developed at the World Radiation Center Physikalish-Meteorologisches Observatorium in Davos, Switzerland (WRC/PMOD) between 2001and 2004(±0.01, Wehrli, 2000)), and finally with CIMEL sun photometer within AERONET (Holben et al., 1998) since 2005 onwards (±0.01,Eck et al., 1999).In this study, we have considered daily AOD values at 500 nm.For 1992 and 1993 there are no AOD data, so we considered a constant AOD value of 0.02 that represents ∼ 80 % of the days at IZO.Estimated GSR using an AOD of 0.01 and 0.03 shows differences with estimated GSR using AOD = 0.02 of 0.32 % and −0.18 %, respectively.

Sunshine duration data
SD is the time period that the ground surface is irradiated by direct solar radiation (i.e.sunlight reaching the Earth's surface directly from the sun).In 1982, WMO defined it as the period during which direct solar irradiance exceeds a threshold value of 120 W m −2 (WMO, 1982).This value accounts for the expected direct solar irradiance shortly after (before) sunrise (sunset) in cloud-free, but also in low aerosol load conditions (WMO, 1984).At IZO, the SD observations started in 1933 with a Campbell-Stokes sunshine recorder (onwards, CS).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 the burn when the card is removed from the instrument at the end of the day (Painter, 1981;WMO, 1996;Sanchez-Lorenzo et al., 2013).During recent years, this traditional sunshine recorder has been replaced by electronic devices.At IZO, the CS instrument was replaced by a Kipp and Zonen Sunshine Duration Sensor (CSD) in 2001, which still operates.
In order to document the precision of the IZO SD measurements as observed by the CS, we have compared these measurements to those obtained from DSR simulated with LibRadtran when exceeding a threshold value of 120 W m −2 since DSR measurements are not available during the CS data series period (see Sect. 2.2 for details about the simulations).We have considered all the cloud-free days, selected by using the method of Long and Ackerman (2000), with low aerosol content because they ensure very stable atmospheric conditions at IZO, mainly observed from October to February and between May and June (Rodríguez et al., 2011, and references therein).Notice that the 1997/1999 period has been selected throughout this work as testing period (see Sect. 4.1).
The comparison shows a good agreement (see Fig. 1), obtaining a correlation coefficient and a RMSE of 0.95 and 0.52 h (4.8 %), respectively.Nonetheless, we observe that CS records systematically overestimate the sunshine hours by 3.1 %, showing a seasonal dependence: the CS record tends to overestimate SD by 2.4 % from October to February, and 5.6 % in May and June.This seasonal variation, also found by other authors (e.g.Kerr and Tabony, 2004;Hinssen and Knap, 2007), may partly be attributed to the different response of the CS recorder to atmospheric conditions in R. D. García et al.: Reconstruction of global solar radiation time series at Iza ña station f 0.09 MJm −2 (0.31 %) and 0.04 MJm −2 comparing with solar observations, the mean ations-observations) is -0.30±0.24MJm −2 ) for GSR and -0.16±0.34MJm −2 (for DSR. ine duration data me period that the ground surface is irradiect solar radiation (i.e., sunlight reaching the ace directly from the sun).In 1982, WMO dee period during which direct solar irradiance hreshold value of 120W m −2 (WMO, 1982).accounts for the expected direct solar irradiy after (before) sunrise (sunset) in cloud-free, ow aerosol load conditions (WMO, 1984).the SD observations started in 1933 with -Stokes sunshine recorder (onwards, CS).The the direct solar beam through a glass sphere, ncentrically in a section of a spherical bowl, on located begin the sphere.The card is provided indication, which makes it possible to determine the length of the burn when the card is removed strument at the end of the day (Painter, 1981;;Sanchez-Lorenzo et al., 2013).During the last aditional sunshine recorder has been replaced by evices.At IZO, the CS instrument was replaced nd Zonen Sunshine Duration Sensor (CSD) in operates until now.o document the precision of the IZO SD meas observed by the CS, we have compared these ts to those obtained from DSR simulated with when exceeding a threshold value of 120 W m −2 easurements are not available during the CS period (see Sect. 2.2 for details about the sime have considered all the cloud-free days, seing the method of Long and Ackerman (2000), rosol content because they ensure very staeric conditions at IZO, are mainly observed whether the ambient air is humid (typically in winter) or 315 dry (typically in summer) (Wood et al., 2003) as well as the burning of the card strip is not well defined at sunrise and sunset, leading to differences through the year.Also, the uncertainties introduced by the model, about 1 %, should be taken into account (see Sect. 2.2).

3 Recalibration of measured GSR
The GSR measurements between January 1992 and July 2005 were not performed following strict quality controls.Therefore, a recalibration of the corresponding data is mandatory in order to be reliably used in the reconstruction 325 of the IZO GSR time series.The recalibration was done by using the LibRadtran, which has demonstrated to be a very useful quality control tool (García et al., 2014).To do so, we calculated the new sensitivity of our pyranometers (here- winter and summer months.For example, the card strip reacts in a different manner depending on whether the ambient air is humid (typically in winter) or dry (typically in summer) (Wood et al., 2003) and also the burning of the card strip is not well defined at sunrise and sunset, leading to differences through the year.Further, the uncertainties introduced by the model, about 1 %, should be taken into account (see Sect. 2.2).

Recalibration of measured GSR
The GSR measurements between January 1992 and July 2005 were not performed following strict quality controls.Therefore, a recalibration of the corresponding data is mandatory in order for it to be reliably used in the reconstruction of the IZO GSR time series.The recalibration was done by using the LibRadtran, which has demonstrated to be a very useful quality control tool (García et al., 2014).To do so, we calculated the new sensitivity of our pyranometers (hereafter re-evaluated sensitivity,  interval of 5 min from the solar noon (WMO, 1996).These GSR estimations were done for all the cloud-free days with low aerosol content in the months from October to February and between May and June to assure very stable atmospheric conditions.This approach was tested with BSRN measurements between 2009 and 2012, obtaining a difference < 3 % between the original sensitivity given by the manufacturer and the re-evaluated one (García et al., 2014).
The recalibrated daily GSR measured at IZO and OT between 1992 and 2005 is illustrated in Fig. 2. The relative difference between re-evaluated and original calibration is, on average, < 7 %, except from April 1997 to August 1998 (Fig. 2a).In this period the re-evaluated calibration decreases by 15 % compared to the original calibration, with the consequent increase of daily GSR (see Fig. 2b).Fig. 2b and c show the time series of daily GSR from 1992 to 2005 at IZO (red squares represent original GSR and blue squares represent re-evaluated GSR) and OT (red squares original GSR and green squares re-evaluated GSR), respectively, and the whole recalibrated time series is plotted in Fig. 2d.Note that the global GSR measured from SCP instrument was extended to 2800 nm with LibRadtan to be comparable to the wavelength range of the other pyranometers.

Estimation of GSR from sunshine duration
Several types of regression models have been proposed for estimating GSR on a horizontal surface at the Earth's surface from SD records.One of the most extended and used approaches was developed by Ångström (Angstrom, 1924(Angstrom, , 1956) ) and later modified by Prescott and Rietveld (Prescott, 1940).This model provides GSR from SD by using the following equation: where H o is the extraterrestrial solar radiation on a horizontal surface (MJ m −2 day −1 ), n and N d (Eq.2) are the number of hours measured by the SD recorder and the maximum daily SD, respectively, and a and b are empirically determined regression constants: where φ is the geographic latitude and δ is the solar declination (Eq.5) The value of H o is calculated as where I sc is the solar constant, E o is the eccentricity correction factor of the Earth's orbit (Eq.4) and ω s is the sunrise hour angle (Spencer, 1971).Note that we have used monthly values of the solar constant instead of a constant value: where d n is the day number of the Julian day of the year, starting from the first of January.
The SD records and, thus, the fraction of clear sky (FCS) defined here by Eq. ( 7), depend on solar direct irradiance, cloudiness (amount, type and thickness), PWV and atmospheric aerosols (mainly mineral dust particles at IZO, Rodríguez et al., 2011;O. E. García et al., 2012): They also depend, to a lesser extent, on meteorological variables as temperature and humidity, although this dependence is very small and purely instrumental.All of these factors account for the stratification found in Fig. 3a, where five regions of FCS values (in intervals of 20 %) can clearly be distinguished, with a very low overlapping among them.Similar stratification is observed in the measured GSR time series (Fig. 3b).Therefore, the subsequent estimation of GSR from SD records, by using Ångström-Prescott's Eq. (1), will be performed considering the dependence on FCS.In the next sections we will evaluate the accuracy of our approach by comparing the GSR estimations from SD records with GSR measurements between 1992 and 2000 when simultaneous GSR and SD measurements are available.

Validation of GSR estimations
In order to validate our GSR estimations, we have split the SD series (1992/2000) into two periods: one for calculating the Ångström-Prescott coefficients ("calculation period"), and another for testing them ("testing period").Once the coefficients a and b are evaluated in the calculation period, we compute daily GSR estimates using the Eq. ( 1) for the testing period.
As calculation period we consider the period 1992/1996 with high-quality GSR measurements to assure very reliable coefficients, (when minor post-corrections were applied, see Fig. 2), and 1997/1999 was used as testing period.In accordance with Fig. 3, the Ångström-Prescott coefficients have been computed by grouping FCS values into 20 % intervals.Table 3 lists the obtained coefficients.
Theoretically, the expected uncertainty on the GSR estimations can be determined by doing an error propagation dependence is very small and purely instrumental.All of these factors account for the stratification found in Fig. 3a, where five regions of FCS values (in intervals of 405 20 %) can clearly be distinguished, with a very low overlapping among them.Similar stratification is observed in the measured GSR time series (Fig. 3b).Therefore, the subsequent estimation of GSR from SD records, by using Ångström-Prescott's Eq. (1), will be performed considering 410 the dependence on FCS.In the next sections we will evaluate the accuracy of our approach by comparing the GSR estimations from SD records with GSR measurements between 1992 and 2000 when simultaneous GSR and SD measurements are available.

Validation of GSR estimations
In order to validate our GSR estimations, we have split the SD series (1992/2000) into two periods: one for calculating the Ångström-Prescott's coefficients ("calculation period"), and another for testing them ("testing period").Once 420 the coefficients a and b are evaluated in the calculation period, we computed daily GSR estimates using the Eq. ( 1) for the testing period.
As calculation period we consider the period 1992/1996 with high-quality GSR measurements to assure very reliable 425 coefficients, (when minor post-corrections were applied, see Fig. 2), and 1997/1999 was used as testing period.In accordance with Fig. 3, the Ångström-Prescott's coefficients have been computed by grouping FCS values into 20 % intervals.Table 3 lists the obtained coefficients.

430
Theoretically, the expected uncertainty on the GSR estimations can be determinated by doing an error propagation on Eq. (1).Thus, assuming the errors given in Table 3 for the Ångström-Prescott's coefficients and the errors given in  on Eq. ( 1).Thus, assuming the errors given in Table 3 for the Ångström-Prescott coefficients and the errors given in Sect.2.3 for the SD measurements, we obtain that the GSR estimations can be provided with an error of 4.0 % for the corresponding testing period.
This theoretical quality estimation has been complemented with a detailed experimental intercomparison summarised in Fig. 4 and Table 4.The latter shows the statistics for the bias between GSR estimations and measurements (estimations − measurements) at different timescales (daily, monthly, seasonal and annual) and on the FCS intervals for the testing period, whereas the former displays the annual cycle of the bias and its dependence on the FCS intervals.The intra-annual bias reveals that the GSR estimations (Fig. 4a) are more accurate in summer (RMSE of 3.1 %) than in winter (RMSE of 9.3 %), when the meteorological conditions are more variable.This is in line with the observed dependence on the FCS values (Fig. 4b): the lower the FCS values (i.e.high presence of clouds or aerosols) that are observed, the larger the systematic bias and scattering found.This leads to a decrease of correlation between estimations and measurements (Table 4).
The overall systematic bias (median bias) is estimated to be < 0.30 MJ m −2 , while the precision, given by the RMSE, reaches 2.16 MJ m −2 (9.2 %) (see Table 4).However, when considering FCS values > 40 % (∼ 90 % of the days), the RMSE values decrease to 4.5 %.These values agree with our theoretical error estimation and are comparable to previous studies.Several authors reported RMSE from the comparison of the observed GSR and the estimated GSR using the Ångström-Prescott method of 1.26 MJ m −2 in the city of Toledo, Spain (Almorox et al., 2005), between 1.49 and 1.65 MJ m −2 in Ankara, Turkey (Yorukoglu and Celik, 2006) and between 1.39 and 3.08 MJ m −2 in 31 sites around China (Che et al., 2005).

Long-term consistency of GSR estimations
In order to reliably use our GSR estimations from SD recorder for solar radiation trends studies, it is indispensable to document their homogeneity and long-term stability.
To do so, we examine possible drifts and discontinuities in the times series of the differences between the GSR estimations and measurements in the common period (1992/2000).We defined a drift as the linear trend of the monthly median bias (estimations − measurements), while the change-points (changes in the median of the bias time series) are analysed by using a robust rank order change-point test (Lanzante, 1996).Note that we have applied the Ångström-Prescott coefficients obtained during the period 1992-2000 (see Table 5).
The straightforward comparison between the monthly median anomalies reveals a rather consistent agreement (correlation coefficient of 0.66, Fig. 5a), although a systematic change point was detected in the median bias time series in June 1995 at 99 % confidence level (Fig. 5b).This change point may be caused by a change of the glass sphere of the CS recorder in that period.Nevertheless, we observe that there are no significant drifts in the bias time series before (+0.012 ± 0.019 MJ m −2 year −1 linear trend) and after (−0.006 ± 0.013 MJ m −2 year −1 linear trend) this systematic change point at 99 % of confidence level.These findings indicate that the GSR estimates are consistent over time with GSR measurements and that they are valid for reconstructing the global GSR time series and trends studies.ependence is very small and purely instrumental.All f these factors account for the stratification found in ig.3a, where five regions of FCS values (in intervals of 0 %) can clearly be distinguished, with a very low overapping among them.Similar stratification is observed in he measured GSR time series (Fig. 3b).Therefore, the ubsequent estimation of GSR from SD records, by using ngström-Prescott's Eq. (1), will be performed considering he dependence on FCS.In the next sections we will evaluate he accuracy of our approach by comparing the GSR estiations from SD records with GSR measurements between 992 and 2000 when simultaneous GSR and SD measureents are available.
.1 Validation of GSR estimations n order to validate our GSR estimations, we have split he SD series (1992/2000) into two periods: one for calcuating the Ångström-Prescott's coefficients ("calculation peiod"), and another for testing them ("testing period").Once he coefficients a and b are evaluated in the calculation peiod, we computed daily GSR estimates using the Eq. ( 1) for he testing period.
As calculation period we consider the period 1992/1996 ith high-quality GSR measurements to assure very reliable oefficients, (when minor post-corrections were applied, ee Fig. 2), and 1997/1999 was used as testing period.In ccordance with Fig. 3, the Ångström-Prescott's coefficients ave been computed by grouping FCS values into 20 % inervals.Table 3 lists the obtained coefficients.
Theoretically, the expected uncertainty on the GSR estiations can be determinated by doing an error propagation n Eq. (1).Thus, assuming the errors given in Table 3 for he Ångström-Prescott's coefficients and the errors given in Sect.2.3 for the SD measurements, we obtain that the GSR 435 estimations can be provided with an error of 4.0 % for the corresponding testing period.This theoretical quality estimation has been complemented with a detailed experimental intercomparison summarized in Fig. 4 and Table 4.The latter shows the statis-440 tics for the bias between GSR estimations and measurements (estimations − measurements) at different time scales (daily, monthly, seasonal and annual) and on the FCS intervals for the testing period, whereas the former displays the annual cycle of the bias and its dependence on 445 the FCS intervals.The intra-annual bias reveals that the GSR estimations (Fig. 4a) are more accurate in summer (RMSE of 3.1 %) than in winter (RMSE of 9.3 %), when the meteorological conditions are more variable.This is in line with the    (correlation coefficient of 0.66, Fig. 5a), although a systematic change point was detected in the median bias time series in June 1995 at 99% confidence level (Fig. 5b).(Lanzante, 1996) on the annual mean anomalies time series, we confirm three change-points in 1953, 1968 and 2000 at 99 % confidence level.The change-point signal-to-noise ratios, which quantify the magnitude or relative importance of each discontinuity, are 2.3, 0.8 and 1.5, respectively, indicating that the breaks in the 1950s and 2000s may be considered as the principal transition dates and 1968 as a secondary change point.These two principal discontinuities define three periods: 1. Early brightening: from the 1930s to the early 1950s, the few historic data available suggest an increase of the GSR in the first part of the 20th century, known as early brightening (De Bruin et al., 1995;Gilgen et al., 1998;Ohmura, 2006), with a peak around the late 1940s/early 1950s.This phenomenon is also confirmed by the IZO anomalies time series considering all-sky conditions, and especially visible for mostly cloud-free situations (FCS > 40 %, Fig. 7b).However, we observe a delay of between 5 and 10 years for the transition from the early brightening to the dimming period.This delay is generally observed throughout the whole anomalies time series and may be partly attributed to the free troposphere conditions of IZO.Note that between 1933 and 1953 the cloudy and cloud-free anomalies time series show opposite trends, with an anti-correlation of ∼ 25 % (Fig. 7).
2. Dimming: from the 1950s to the end of the 1990s a gradual decrease of GSR is observed, which is in accordance with the widespread period of reduced solar radiation at a global scale, extensively reported by the literature and known as dimming (Stjern et al., 2009;Gilgen et al., 1998;Stanhill and Cohen, 2001;Wild et al., 2005;Ohmura, 2006;Wild, 2009).and Fernandina Island (1968), Chinchón (1982) and Pinatubo (1991).
The dimming/brightening phenomena have consistently been detected both under cloudy and cloud-free conditions, as shown in Fig. 7.This suggests that the dimming and brightening phenomena might be associated with, among others, changes on atmospheric aerosols concentrations due to anthropogenic activities (Wild, 2009, and references herein).Note that cloudy and cloud-free anomalies time series show an acceptable agreement from the 1960s onwards, with a positive and significant correlation of ∼ 50 %.

Summary and conclusions
The GSR times series between 1933 and 2013 has been successfully reconstructed combining GSR estimates from SD measurements (1933/1991) using the Ångström-Prescott method and GSR observations performed by different pyranometers (1992/2013) at IZO.
The quality of the GSR and SD databases have been assessed.Since 2005, the global GSR measurements at IZO are taken following strict quality controls in the framework of the NRC and BSRN networks.Nonetheless, before 2005, recalibration of the GSR measurements was needed.This recalibration was made with the LibRadtran, obtaining differences of < 7 % between the original and the re-evaluated calibration factors.
The SD measurements taken by Campbell-Stokes recorders (CS) between 1993 and 2000 have been validated against the LibRadtran model DSR simulations for cloudfree days.A good agreement is found with a RMSE of 0.52 h (4.8 %).Propagating these errors to GSR estimations from SD records, an expected precision of 4.0 % has been found.By comparing with GSR, we obtain a precision, given by the RMSE, of 2.16 MJ m −2 (9.2 %).However, when considering FCS values > 40 %, the RMSE values drop to 4.5 %.The consistency over time series of the GSR estimations has been documented by comparing these with GSR measurements on a 9-year period (1992/2000).
The resulting annual GSR time series confirms a period of early brightening from the 1930s to the early 1950s, a period of dimming from the 1950s to the ending of the 1990s, followed by a period of brightening in the most recent decades.All these findings demonstrate the consistency of the IZO GSR time series presented in this work, which may be a reference for solar radiation studies in the subtropical North Atlantic region.Future works will analyse in depth the longterm trends and their interplay with interannual variations of the solar constant, cloud cover and the atmospheric aerosols content.The joint analysis with dust AOD is critical in our region since Saharan dust intrusions, which undergo interannual and decadal variations, modulates AOD and hence solar radiation.

Fig. 1 .
Fig. 1.Scatterplot of the measured SD from CS instrument versus SD obtained from LibRadtran model when DSR exceeds a threshold value of 120 Wm −2 .The root mean square error (RMSE), and least-square fit parameters are shown in the legend.

Figure 1 .
Figure 1.Scatter plot of the measured SD from CS instrument versus SD obtained from LibRadtran model when DSR exceeds a threshold value of 120 W m −2 .The root mean square error (RMSE), and least-square fit parameters are shown in the legend.

Figure 2 .
Figure 2. (a) Evolution of the ratio between the re-evaluated calibration and the original calibration for the different pyranometers installed at IZO (blue squares) and at OT (green squares).(b), (c) and (d): time series of the GSR (MJ m −2 ) from 1992 to 2005 at IZO, at OT and the whole recalibrated time series, respectively (red squares represent original GSR, blue squares represent re-evaluated GSR at IZO and green squares represent re-evaluated GSR at OT).

Figure 3 .
Figure 3.Time series of (a) the sunshine duration (hours).(b) GSR from 1992 to 2000 at IZO.The colour scale indicates the fraction of clear-sky values (FCS, %).
Time series of (a) the sunshine duration (hours).(b) GSR rom 1992 to 2000 at IZ0.The color scale indicates the fraction of lear sky values (FCS,%).

Fig. 4 .
Fig. 4. Box plot of the (a) annual cycle of the bias (MJm −2 ) and (b) bias versus fraction of clear sky values (FCS) from 1997 to 1999 at IZO. Lower and upper boundaries for each box are the 25 and 75 percentiles, the solid line is the median value, the crosses indicate values out of the 1.5 fold box area (outliers).

Figure 4 .
Figure 4. Box plot of the (a) annual cycle of the bias (MJ m −2 ) and (b) bias versus fraction of clear sky values (FCS) from 1997 to 1999 at IZO. Lower and upper boundaries for each box are the 25 and 75 percentiles, the solid line is the median value and the crosses indicate values out of the 1.5-fold box area (outliers).

Fig. 5 .
Fig. 5. Times series of monthly median of (a) the deseasonalized anomalies of GSR estimations (red line) and measurements (black line) and (b) monthly median bias between GSR estimations and measurements (MJm −2 ) from 1992 to 2000 at IZO.The error bars indicate ±1 SEM (standard error of the monthly median) and the black arrow indicates the change point date.

Figure 5 .Fig. 6 .Fig. 7 .
Figure 5. Times series of monthly median of (a) the deseasonalized anomalies of GSR estimations (red line) and measurements (black line) and (b) monthly median bias between GSR estimations and measurements (MJ m −2 ) from 1992 to 2000 at IZO.The error bars indicate ±1 SEM (standard error of the monthly median) and the black arrow indicates the change-point date.

Figure 6 .
Figure 6.(a) Time series of the GSR (MJ m −2 ) monthly means.The black dots represent the GSR obtained from SD data and the grey triangles represent the GSR measured between 1992 and 2000.(b) Annual means of the GSR anomalies.The arrows indicate the change-point dates.The error bars indicate ±1 SEM (standard error of the mean).Five-year moving average is shown in red.

Fig. 6 .
Fig. 6.Time series of the GSR (MJm ) (a) monthly means.The black dots represent the GSR obtained from SD data and the gray triangles represent the GSR measured between 1992 and 2000 (b) annual means of the GSR anomalies.The arrows indicate the change-point dates.The error bars indicate ±1 SEM (standard error of the mean).Five-yr moving average is shown in red.

Table 1 .
Pyranometers installed between 1992 and 2013 at IZO and OT.

Table 3 .
Coefficients a and b between 1992 and 1996 for five FCS intervals.

Table 4 .
Statistics for the differences between simulations and measurements at IZO (in MJ m −2 ) between 1997 and 1999 at different timescales (daily, monthly, seasonal and annual) and on the fraction of clear sky (FCS).(RMSE: root mean square error; R: correlation coefficient).The statistics for the relative bias are in parentheses (%).(DJF: December-January-February; MAM: March-April-May; JJA: June-July-August and SON: September-October-November).
(Almorox et al., 2005)nstruction of global solar radiation time series at Iza ña station 7 le 4. Statistics for the differences between simulations and meaements at IZ0 (in MJm −2 ) between 1997 and 1999 at different e scales (daily, monthly, seasonal and annual) and on the fracof clear sky (FCS).(RMSE:RootMeanSquareError;R:correon coefficient.The statistics for the relative bias are in brackets .(DJF:December-January-February;MAM:March-April-May;:June-July-AugustandSON: September-October-November)., while the precision, given by the SE, reaches 2.16 MJ m −2 (9.2 %) (see Table4).Howr, when considering FCS values > 40 % (∼ 90 % of the s), the RMSE values down to 4.5 %.These values agree h our theoretical error estimation and are comparable to vious studies.Several authors reported RMSE from the parison of the observed GSR and the estimated GSR ng Ångström-Prescott method of 1.26 MJ m −2 in the of Toledo, Spain(Almorox et al., 2005), between 1.49 1.65 MJ m −2 in Ankara, Turkey(Yorukoglu and Celik,

Table 5 .
Coefficients a and b between 1992 and 2000 at IZO as a function of the fraction of clear sky values (FCS).SEM is the standard error of the mean.

Table 5 .
Coefficients a and b between 1992 and 2000 at IZO as a function of the fraction of clear-sky values (FCS).SEM is the standard error of the mean.