Global radiative flux and cloudiness variability for the period 1959-2010 in Belgium: A comparison between reanalyses and the regional climate model MAR

The use of regional climate models (RCMs) can partly reduce the biases in global radiative flux (Eg↓) that are found in reanalysis products and global models, as they allow for a finer spatial resolution and a finer parametrisation of surface and atmospheric processes. In this study, we assess the ability of the MAR («Modèle Atmosphérique Régional») RCM to reproduce observed changes in Eg↓, and we investigate the added value of MAR with respect to reanalyses. Simulations were performed at a horizontal resolution of 5 km for the period 1959–2010 by forcing MAR with different reanalysis products: ERA40/ERA-interim, NCEP/NCAR-v1, ERA-20C, and 20CRV2C. Measurements of Eg↓ from the Global Energy Balance Archive (GEBA) and from the Royal Meteorological Institute of Belgium (RMIB), as well as cloud cover observations from Belgocontrol and RMIB, were used for the evaluation of the MAR model and the forcing reanalyses. Results show that MAR enables largely reducing the mean biases that are present in the reanalyses. The trend analysis shows that only MAR forced by ERA40/ERA-interim shows historical trends, which is probably because the ERA40/ERA-interim has a better horizontal resolution and assimilates more observations than the other reanalyses that are used in this study. The results suggest that the solar brightening observed since the 1980s in Belgium has mainly been due to decreasing cloud cover.


Introduction
The amount of solar radiation reaching the Earth's surface or global radiative flux (hereafter referred to as E g↓ ) governs a wide range of processes that support life on Earth. It generates changes in the water cycle and carbon cycle, and causes droughts or heat waves [1]. Furthermore, renewable energies, such as solar power, are a growing part of the energy sector. Thus, the supply system could be more vulnerable to E g↓ changes [2]. Given the impacts on our environment and our electricity production, an accurate insight into changes in the E g↓ and the causes of such changes is crucial.
Measurements provide the best estimate of E g↓ variability where they are taken. Worldwide monitoring of E g↓ from ground-based stations began in the late 1950s. These ground-based measurements suggest that the period 1950-1980 was marked by a decrease in E g↓ that is better known as "global dimming" [3][4][5][6]. This dimming was followed by a partial recovery ("brightening")

The MAR Model
The RCM used in this study is the MAR model coupled to the one-dimensional (1D) surface transfer scheme SISVAT (Soil-Ice-Snow-Vegetation-Atmosphere-Transfer). A more complete description of the atmospheric module of MAR can be found in Gallée and Schayes [35], Gallée [36], and Gallée et al. [37], while the surface transfer scheme is detailed in De Ridder and Gallée [38]. MAR is a hydrostatic primitive equation model in which convection is parametrised according to Bechtold et al. [39]. Cloud microphysical parametrisations are based on the studies of Kessler [40], Lin et al. [41], Meyers et al. [42], and Levkov et al. [43]. The radiative transfer through the atmosphere uses the radiation scheme from ERA-40 based on Morcrette [44]. It uses time-varying aerosols from the monthly climatology of tropospheric aerosols (black carbon, organic, sea salt, sand, and sulphate aerosols) by Tegen et al. [45] and the GISS (Goddard Institute for Space Studies) volcanic aerosols, while the cloud microphysics module uses fixed aerosol concentration [42]. Therefore, we can assume that only the aerosol-radiation interactions and the cloud-radiation interactions are taken into account in MAR.
Although initially designed for polar regions [35,46], MARv3.6 was successfully calibrated to the climate of Belgium [47]. However, by comparing the MAR results to observations from 20 weather stations spread over the Belgian territory, Wyard et al. [45] showed that MARv3.6 is too warm (+1.5 • C) and too dry (−10%) in summer, while it is too cold (−0.5 • C) and too wet (+20%) in winter, as a result of an underestimation of the cloud cover in summer and overestimation in winter. Despite the aforementioned biases in the simulated cloud cover, Hofer et al. [48] concluded that "cloud cover from MAR shows good agreement in the spatial distribution and amplitude of the changes observed over Greenland by both satellite platforms [the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor and (ii) the Advanced Very High Resolution Radiometer (AVHRR)]".
In the present study, we use version 3.8 of MAR. Compared with MARv3.6, which was used in Wyard et al. [47], MARv3.8 reduces temperature biases by 20% by correcting in large parts the underestimation of cloud cover by delaying the onset of precipitation, and hence increasing cloud cover. In addition, the convective scheme is called twice as often in MARv3.8 than in MARv3.6, so that convective clouds reside longer in the atmosphere before precipitating [49]. Finally, in MARv3.8, the vegetation seasonality is also better taken into account by using daily leaf area index (LAI) and green leaf fraction (GLF) values from the MERRA-2 reanalysis [50] compared with MARv3.6, which uses monthly values.
The integration domain (120 × 110 pixels, centred on (50.4 • N; 4.3 • W)), the space grid (5 km), and the number of vertical levels (24 levels between 3-18 km above the surface) are the same as those Atmosphere 2018, 9,262 4 of 23 used in Wyard et al. (2017) [47]. With the aim of evaluating the sensitivity of the MAR results on the lateral boundary conditions, we have forced MAR using four different reanalyses.

Reanalyses
The lateral boundary conditions [temperature, wind, humidity, and pressure at each vertical level of MAR, as well as the sea surface temperature (SST)] were provided every 6 h by nesting MAR in four reanalysis products: (1) ERA-Interim (horizontal resolution 0.75 • ; 60 vertical levels from the surface to 0.01 hPa) from the European Centre for Medium-Range Weather Forecasts (ECMWF), which is a third generation reanalysis assimilating most of the available in situ and satellite data [23]. Wyard et al. [47] showed that MAR forced by this reanalysis provides the best daily mean temperature and precipitation amount estimates over Belgium. ERA-Interim data is available from 1979, but we extended this time series by using ERA40 for the period 1958-1978 [26]. ERA40 (horizontal resolution 1.125 • ; 60 vertical levels) is a second-generation reanalysis, and was the first to assimilate satellite data. The MAR simulation using ERA-Interim/ERA40 as forcing is hereafter referred to as MAR-ERA. (2) NCEP/NCAR-v1 (horizontal resolution 2.5 • ; 28 vertical levels from the surface to 0.3 hPa) from the NCEP-NCAR [24]. This first-generation reanalysis assimilates aircraft, rawinsonde, ship, and station data, as well as satellite retrievals. NCEP/NCAR-v1 data are available from 1948. The MAR simulation using NCEP/NCAR-v1 as forcing is hereafter referred to as MAR-NCEP1. (3) ERA-20C (horizontal resolution 1.125 • ; 91 vertical levels from the surface to 0.01 hPa) is the first atmospheric reanalysis by ECMWF covering the entire 20th century [25]. This third-generation reanalysis only assimilates surface observations such as pressure and marine winds. ERA-20C is available for 1900-2010. The MAR simulation using this forcing is hereafter referred to as MAR-ERA-20C. (4) 20CRV2C (horizontal resolution 2 • ; 28 vertical levels from the surface to 0.3 hPa), the 20th century reanalysis of the National Oceanic and Atmospheric Administration (NOAA) [22]. It is based on an ensemble mean of 56 members assimilating only surface pressure, monthly sea surface temperature, and sea ice cover. 20CRV2C is available for 1900-2010. The MAR simulation using it as forcing is hereafter referred to as MAR-20CRV2C.
MAR runs at 5 km were directly forced by these reanalyses at the boundaries. Finally, the simulations were performed for the period 1959-2010, which is the longest period common to all of the reanalyses considered in this study.

Evaluation Datasets
Measurements of E g↓ and TCC from ground stations over Belgium are used as reference datasets for the evaluation of MAR. E g↓ data were supplied by GEBA and by the Royal Meteorological Institute of Belgium (RMIB). The GEBA E g↓ dataset consists of monthly mean E g↓ , which was measured by well-calibrated instruments and underwent rigorous quality checks [3,51]. The relative random error of GEBA data was found to be approximately 2% of the annual mean values [3]. Four time series consisting of at least 30 years are available for Belgium in the GEBA data set, as listed in Table 1. The RMIB has a long-term experience with ground-based measurements of E g↓ in Belgium with uninterrupted measurements at Uccle since 1951 [52]. Uccle is one of the 22 Regional Radiation Centres established within the World Meteorological Organisation (WMO) regions. In addition to Uccle, RMIB provided E g↓ data from three other stations (Table 1). These data consist of monthly mean E g↓ calculated from 30 min-average measurements. The time-series from GEBA are the same as those from RMIB, with the difference between them lower than 1%. Since the RMIB time-series are longer, we removed these redundant GEBA time-series from our dataset, so that we finally used E g↓ data from the four RMIB stations shown in Figure 1: Melle, Middlekerk/Oostende, Saint-Hubert, and Uccle. Seasonal and annual mean E g↓ were then computed using the monthly data. However, if more than one month was missing within a year, the annual mean was not computed; otherwise, the mean E g↓ for the year was computed by filling the missing month by the climatic average of this month for the period covered by the time series. The seasonal mean E g↓ values are not computed if there are gaps in the monthly data. Table 1. Available E g↓ and total cloud cover (TCC) time series of more than 30 years for Belgium. The percentage of gaps in the monthly data is also indicated. The time-series that were not selected as reference datasets for the evaluation of the Modèle Atmosphérique Régional (MAR) outputs and the reanalysis products have a grey background. GEBA: Global Energy Balance Archive; RMIB: Royal Meteorological Institute of Belgium.

Station
Coordinates  Figure 1). Monthly, seasonal, and annual mean TCC were calculated from hourly or three-hourly observations carried out by human eye inspection of the sky during the day and at night. Missing values were not taken into account, and a maximum of 5% of missing data is tolerated in our computations. Table 1. Available Eg↓ and total cloud cover (TCC) time series of more than 30 years for Belgium. The percentage of gaps in the monthly data is also indicated. The time-series that were not selected as reference datasets for the evaluation of the Modèle Atmosphérique Régional (MAR) outputs and the reanalysis products have a grey background. GEBA: Global Energy Balance Archive; RMIB: Royal Meteorological Institute of Belgium.

Station
Coordinates   The TCC data for four stations were supplied by the RMIB and the Belgian agency managing the national civil airspace (Belgocontrol) ( Table 1 and Figure 1). Monthly, seasonal, and annual mean TCC were calculated from hourly or three-hourly observations carried out by human eye inspection of the sky during the day and at night. Missing values were not taken into account, and a maximum of 5% of missing data is tolerated in our computations.

Comparison of Annual and Seasonal E g↓ and TCC Observations with the MAR Outputs and Reanalyses
In this section, E g↓ and TCC values computed by the different MAR simulations (MAR-ERA, MAR-NCEP1, MAR-ERA20C, and MAR-20CRV2C) and the corresponding reanalyses forcing MAR (ERA40/ERA-interim, NCEP/NCAR-v1, ERA-20C, and 20CRV2C) are compared to the E g↓ and TCC measurements obtained from various ground stations (Table 1 and Figure 1). As all of the stations show similar results, we only display here the statistics computed from E g↓ and TCC measured at Saint-Hubert, because it has the most complete time-series for both investigated parameters.

Comparison of Annual and Seasonal Eg↓ and TCC Observations with the MAR Outputs and Reanalyses
In    Figure 4d). Furthermore, the comparison suggests that MAR slightly overestimates E g↓ with mean annual biases (MB) ranging from +7.7 W m −2 (MAR-NCEP1) to +9.9 W m −2 (MAR-ERA) (Figure 4b). On a seasonal time-scale, we find the lowest values of R in the summer (June-July-August, JJA) ranging from 0.62 (MAR-20CRV2C) to 0.70 (MAR-ERA), while CRMSE was among the largest, with values above 13.8 W m −2 . The other seasons exhibit values of R of at least 0.72 (MAR-NCEP1 in autumn (September-October-November, SON)). The lowest CRMSE was found in winter (December-January-February, DJF) with values below 3.9 W m −2 (MAR-NCEP1). MAR tends to overestimate E g↓ for all of the seasons except winter. We find the largest MB in the summer, with values ranging from +22.3 W m −2 (MAR-ERA-20C) to +29.2W m −2 (MAR-20CRV2C). In winter, MB is negative, and ranges from −3.  Finally, we observe the same general behaviour for Eg↓ modelled by MAR compared to Eg↓ directly from the driving reanalyses. The values of R between reanalyses and measurements are as high as between MAR and the measurements. MB, for its part, is largely reduced in MAR compared to the reanalyses (Figures 2b and 4b). For instance, MAR-NCEP1 shows an annual MB of +7.7 W m −2 while this value reaches +64.3 W m −2 in NCEP/NCAR-v1, which is in agreement with the findings of Zhang et al. [31]. We also find differences between ERA40 (1962ERA40 ( -1978 and ERA-interim (1979-2010), which are fortunately not observed in the associated MAR simulation (MAR-ERA) ( Figure 2). Indeed, ERA40 systematically underestimates Eg↓ (annual MB = −5.00 W m −2 ), while ERA-Interim overestimates Eg↓ (annual MB = +15.6 W m −2 ) in agreement with Zhang et al. [31]. In addition, the forcing reanalyses and MAR have generally comparable CRMSE values. This is because the interannual variability of MAR is fully driven by the forcings. However, the forcing reanalyses exhibit much larger RMSE than MAR because of their larger biases. Finally, we observe the same general behaviour for E g↓ modelled by MAR compared to E g↓ directly from the driving reanalyses. The values of R between reanalyses and measurements are as high as between MAR and the measurements. MB, for its part, is largely reduced in MAR compared to the reanalyses (Figures 2b and 4b). For instance, MAR-NCEP1 shows an annual MB of +7.7 W m −2 while this value reaches +64.3 W m −2 in NCEP/NCAR-v1, which is in agreement with the findings of Zhang et al. [31]. We also find differences between ERA40 (1962ERA40 ( -1978 and ERA-interim (1979-2010), which are fortunately not observed in the associated MAR simulation (MAR-ERA) ( Figure 2). Indeed, ERA40 systematically underestimates E g↓ (annual MB = −5.00 W m −2 ), while ERA-Interim overestimates E g↓ (annual MB = +15.6 W m −2 ) in agreement with Zhang et al. [31]. In addition, the forcing reanalyses and MAR have generally comparable CRMSE values. This is because the interannual variability of MAR is fully driven by the forcings. However, the forcing reanalyses exhibit much larger RMSE than MAR because of their larger biases. E g↓ , the MAR performances also depend on the season. In the summer, R is the lowest, and does not exceed 0.57 in MAR-ERA. Meanwhile, MB values are negative, and range from −5.2% (MAR-ERA-20C) to −13.9% (MAR-20CRV2C), while the CRMSE of MAR-20CRV2C reaches 9.8% in summer. In winter, the R values are among the largest, reaching 0.72 in MAR-ERA. Values of MB are positive, and range from +3.0% (MAR-ERA-20C) to +8.4% (MAR-NCEP1) with the largest CRMSE in winter reaching 4.8% in MAR-NCEP1. As for E g↓ , this comparison highlights that MAR-NCEP1 shows the worst results. In order to explain this apparent artefact in TCC modelled by MAR-NCEP1, we have compared all three cloud types computed by MAR: low cloud cover (LCC), which is defined as the cloud fraction below 680 hPa, medium cloud cover (MCC), which is defined as the cloud fraction between 680-440 hPa, and high cloud cover (HCC), which is defined as the cloud fraction above 440 hPa ( Figure S7, Supplementary Materials). The results suggest that this artefact likely results from an inaccurate simulation of HCC, and therefore from biases in the forcing fields of MAR in the high troposphere ( Figure S7c). have compared all three cloud types computed by MAR: low cloud cover (LCC), which is defined as the cloud fraction below 680 hPa, medium cloud cover (MCC), which is defined as the cloud fraction between 680-440 hPa, and high cloud cover (HCC), which is defined as the cloud fraction above 440 hPa ( Figure S7, Supplementary Materials). The results suggest that this artefact likely results from an inaccurate simulation of HCC, and therefore from biases in the forcing fields of MAR in the high troposphere ( Figure S7c). The comparison with the statistics computed from the TCC modelled by the reanalyses was uneasy, as TCC data from NCEP/NCAR-v1 and 20CRV2C were not available. However, the TCC modelled by MAR shows comparable R, MB, and CRMSE, but MAR shows larger CRMSE values than the corresponding reanalyses. The inhomogeneity between ERA40 and ERA-Interim between 1978 and 1979 is also present in TCC. Regardless of the season (except in the summer), the MB is The comparison with the statistics computed from the TCC modelled by the reanalyses was uneasy, as TCC data from NCEP/NCAR-v1 and 20CRV2C were not available. However, the TCC modelled by MAR shows comparable R, MB, and CRMSE, but MAR shows larger CRMSE values than the corresponding reanalyses. The inhomogeneity between ERA40 and ERA-Interim between 1978 and 1979 is also present in TCC. Regardless of the season (except in the summer), the MB is systematically positive in ERA40, while it is negative in ERA-Interim at Saint-Hubert (Figures 3b and 5a).
Biases in E g↓ are in agreement with the biases found in TCC. For instance, the underestimation of E g↓ in winter is in agreement with the overestimation of TCC found in the models. It is the opposite in summer. Winter exhibits the smallest CRMSE and the highest R for both variables, while summer exhibits the largest CRMSE and the lowest R. More generally, several studies highlight the difficulties encountered by all of the climate models to represent cloud physics. Markovic et al. [53] showed that ERA40 has a tendency to underestimate very low values of E g↓ during the winter, and suggested that this might be due to an underestimation of the reflectivity of optically thick clouds during winter in ERA40. Träger-Chatterjee et al. [30] showed that cloud physics and E g↓ have not been improved in ERA-Interim compared with ERA40, despite improvements in the water cycle representation, model resolution, and data assimilation system. Zhang et al. [31] showed that reanalyses tend to overestimate E g↓ as a result of underestimation in cloud fraction. Alexandri et al. [54] showed that cloud fractional cover, cloud optical thickness, and aerosol optical depth are the most important factors that explain the biases in E g↓ that are found in the RegCM4 RCM over the European domain. Finally, Wyard et al. [47] showed that MAR underestimates convective clouds with consequences on summer temperature (overestimation) and precipitation (underestimation). Given that MAR simulates smaller E g↓ and TCC values than the reanalyses, which overestimate both E g↓ , and TCC, we assume that the cloud optical depth (COD) is larger in MAR than in the reanalyses. Furthermore, part of the biases in E g↓ will also result from the exclusion of the aerosol-cloud interactions in MAR. Fan et al. [55] showed that inaccuracies in the aerosol-cloud interactions contributed to the biases that were found in the GCMs of CMIP5 (Coupled Model Intercomparison Project Phase 5). Finally, the reliability of cloudiness observations is also questionable, given that synoptic observations are performed by human inspection of the sky [56].

Trends in Annual E g↓ and TCC during the Dimming and the Brightening Periods
A linear trend analysis is performed for 1959-1979 and 1980-2010, corresponding to the dimming and the brightening periods reported in the literature [1]. As in Section 3.1, we use the observations and simulations for Saint-Hubert that are shown in Figures 2 and 3. It should be noted that TCC observations are only available from 1966 for Saint-Hubert. Therefore, trends in TCC and E g↓ measured in Saint-Hubert are computed over 1966-1979 for the dimming period. As done by Wyard et al. [47], trend significance is assessed by using the Snedecor uncertainty range for the 95% confidence interval [57] (Appendix A). The uncertainty range represents the interannual variability of the considered time series.
During the dimming period, we find a non-significant negative trend in E g↓ measurements for Saint-Hubert for the period 1966-1979 valued at −12.9 W m −2 decade −1 ( Table 2). During the brightening period, we find a significant positive trend at Saint-Hubert for the period 1980-2010 valued at +6.2 W m −2 decade −1 . The TCC observed in Saint-Hubert shows no significant trend, although synoptic cloud observations are subject to large uncertainties [56]. This nonetheless suggests that cloudiness changes are not the only factor affecting solar radiation variability. Indeed, the RMIB [58] highlighted a significant decrease in E g↓ at Uccle for the period 1951-1984 valued at −6.3%. They attribute this trend to a sustained increase in anthropogenic aerosols. De Bock et al. [59] and the RMIB [58] reported a significant increase in E g↓ at Uccle valued at +5 ± 2 W m −2 decade −1 for 1985-2014. They also associated it to a non-significant decrease in aerosol optical depth at Uccle for 1985-2014 that is valued at −8 ± 5% per decade, which presents a stabilisation after 2002. This decline in aerosol optical depth results from the enormous efforts to curb emissions and air pollution since the 1980s [60]. The most important effort is probably the 1979 Convention on Long-range Transboundary Air Pollution, which was later extended by eight protocols that identify specific measures to be taken by parties to cut their emissions of air pollutants [60].
Regarding E g↓ modelled by MAR, there are no significant trends for 1966-1979. However, trends are negative during 1959-1979, but only MAR-ERA shows a significant trend valued at −10.4 W m −2 decade −1 for this period ( Table 2). This decreasing E g↓ coincides with an increasing trend in the simulated TCC, which is significant in MAR-ERA and valued at +7.1% decade −1 . On the contrary, for 1980-2010, the trends in E g↓ are positive except in MAR-ERA-20C, and significant in MAR-ERA (+3.8 W m −2 decade −1 ) and MAR-NCEP1 (+7.9 W m −2 decade −1 ) as in observations. This increasing E g↓ coincides with a decreasing trend in simulated TCC, except in MAR-ERA-20C. These trends are significant in MAR-ERA and MAR-NCEP1, and are valued at −1.3% decade −1 and −5.9% decade −1 , respectively. The strong negative trend in MAR-NCEP1 seems to partly result from biases in the forcing fields of MAR in the high troposphere before the 1990s as highlighted in Section 3.1 ( Figure S7c). It should be noted that MAR-ERA-20C and MAR-20CRV2C systematically show non-significant and weaker trends than MAR-ERA, MAR-NCEP1, and observations. This suggests that these reanalysis products covering the whole 20th century (ERA-20C and 20CRV2C) are less reliable, possibly because they assimilate surface observations only. We also computed the trends from E g↓ and TCC modelled by the reanalyses. As ERA40/ERA-Interim is biased by the transition from ERA40 to ERA-Interim between 1978 and 1979, we removed the MB of each part of the time-series (1959-1978 and 1979-2010) to homogenise the times series before computing the trend. The signal of the trends in E g↓ modelled by the reanalyses is of the same order as the one modelled by MAR forced by the corresponding reanalyses (Table 2). This suggests that the trends found in the MAR outputs are partially driven by the forcing reanalyses. However, MAR amplifies the trends in E g↓ compared with the observations and the forcing reanalyses. As highlighted in Section 3.1, the cloud optical depth is likely larger in MAR than in the reanalyses. Therefore, changes in TCC impact E g↓ in MAR more than in the reanalyses.

Trends in Seasonal Mean E g↓ and TCC during the Dimming and Brightening Periods
A linear trend analysis is performed for 1959-1979 and 1980-2010 for each pixel of the MAR domain in order to highlight spatial and seasonal variabilities. Trend significance is again assessed by using the Snedecor uncertainty range for the 95% confidence interval. Only the trends from the MAR-ERA outputs are computed, as they fit the best the observations, which is in agreement with Zhang et al. [31] (Sections 3.1 and 3.2).
Over 1959-1979 (dimming period), the MAR-ERA outputs show an overall significant decrease in E g↓ for all of the seasons (Figure 6a-d). The decrease in E g↓ is especially marked in spring (March-April-May, MAM) and summer, and reaches values between −16 and −18 W m −2 decade −1 in some regions of the country (Figure 6b,c). Chiacchio and Wild [9] and Sanchez-Lorenzo et al. [5] found similar results over a larger European domain. MAR-ERA show a significant increase in TCC in all of the seasons Figure 6e-h). This increase in TCC is the most marked in spring, summer, and autumn, and can reach values between +8 and +10% decade −1 in some parts of the country (Figure 6f-h). The patterns of the trends in TCC and in E g↓ are broadly similar, which confirms that changes in cloud cover drive trends in surface solar radiation. To refine our analysis, we also compute trends for all three MAR cloud types (LCC, MCC, and TCC) defined in Section 3.1. These trends show that the increase in TCC for 1959-1979 is mainly caused by a significant increase in LCC (Figure 7a-e) and HCC (Figure 7k-o). It should be noted that HCC is known to affect longwave radiation, mainly while we are studying shortwave radiation changes [61]. MCC shows a significant increase in spring over , with values ranging from +2 and +10% decade −1 (Figure 7h). Over 1980-2010 (brightening period), measurements at Saint-Hubert showed positive significant trends in spring and summer valued at +10.7 W m −2 decade −1 and +12.0 W m −2 decade −1 , respectively (Table S1, Supplementary Materials). Seasonal trends computed from MAR-ERA outputs also show a significant increase in E g↓ in spring and summer, with values between +4 and +10 W m −2 decade −1 (Figure 8b,c). Regarding TCC, the observations show non-significant trends for all seasons (Table S1), while the values modelled by MAR-ERA show a significant decrease in TCC, especially in summer (Figure 8g) reaching −2 to −4% decade −1 over the entire country. Seasonal trends in cloud types computed in MAR-ERA show a significant decrease in LCC for 1980-2010, especially in spring, with values ranging between −2 and −6% decade −1 , depending on the location (Figure 9c). Our results are in agreement with Eastman and Warren [62], who also found a decrease in observed cloudiness in Western Europe for 1971-2009. Previous studies also stated that the brightening results from the combination of both decreasing aerosol emissions and decreasing LCC and MCC [59]. However, after the 1990s, the decrease in cloudiness might have become the dominant factor, according to the results of Mateos et al. [63], who drew similar conclusions over Spain, and the results of Sanchez-Lorenzo et al. [64] and Pfeifroth et al. [65] over Europe.

Origin of Cloudiness Changes
The analysis of the linear trends (Sections 3.2 and 3.3) has highlighted that E g↓ changes are driven by changes in low and medium cloud cover. In this section, we investigated the origin of such changes. As output from MAR is dependent on the reanalysis product that is used to force the model (Section 3.2), we analysed the forcing variables from all four reanalyses at two pressure levels (700 hPa and 500 hPa), which are the characteristic levels of low and medium cloud cover. We especially paid attention to the air temperature (TA), specific humidity (HUS), and relative humidity (HUR). We used the data from the pixel containing Saint-Hubert ( Figure S8, Supplementary Materials), and computed the linear trends and their significance for the periods 1959-1979 and 1980-2010 (Table 3). Table 3. Annual trends and their Snedecor uncertainty range for the 95% confidence interval computed for 1959-1979 and 1980-2010. The values in bold indicate that the trends are statistically significant at the 95% confidence level, namely that they are larger than the associated uncertainty range. For the period 1959-1979 (dimming period), both ERA40/ERA-Interim and NCEP/NCAR-v1 show a significant increase in HUR at 500 hPa valued at +5.4% decade −1 and +2.8% decade −1 , respectively (Table 3). ERA40/ERA-Interim also shows a significant increase in HUR at 700 hPa that was valued at +4.1% decade −1 . This coincides with a positive but non-significant trend in TA at both pressure levels, but also with a negative but non-significant increase in HUS. Only ERA40/ERA-interim shows a significant trend in HUS at 500 hPa valued at +0.05 g kg −1 decade −1 . Both ERA20C and 20CRV2C simulate changes in HUR and TA that are similar to ERA40/ERA-Interim and NCEP/NCAR-v1, but these changes are weaker. In addition, ERA20C and 20CRV2C exhibit positive (but non-significant) trends in HUS at 700 hPa and 500 hPa, while these trends are negative in ERA40/ERA-Interim and NCEP/NCAR-v1 for 1959-1979. Given these trends, increasing LCC and MCC during the dimming period could result from decreasing TA and increasing HUS generating an increase in HUR, which could have therefore enhanced cloud development. The initial changes in TA and HUS could be explained by changes in the atmospheric circulation. In addition, we know that the dimming period was characterised by a high polluted atmosphere [60]. As aerosols suspended in the troposphere could have served as condensation nuclei (aerosol-cloud interactions), they could be responsible for the increase in LCC and MCC during the dimming period by enhancing cloud development. Then, the increase in cloud cover could generate less E g↓ and therefore decreasing TA, which could generate an increase/decrease in HUR.

ERA40/ERA-Interim
For the period 1980-2010 (brightening period), both ERA40/ERA-Interim and NCEP/NCAR-v1 show a significant decrease in HUR at both pressure levels ranging from −1.3% decade −1 (ERA40/ERA-interim) to −1.5% decade −1 (NCEP/NCAR-v1) at 700 hPa, and from −1.1% decade −1 (ERA40/ERA-interim) to −1.5% decade −1 (NCEP/NCAR-v1) at 500 hPa. This coincides with a positive but non-significant trend in TA at both pressure levels, but also with a negative but non-significant increase in HUS. ERA20C and 20CRV2C both exhibit weak non-significant decreasing trends in HUR. These trends coincide with a significant increase in TA at both pressure levels, ranging from +0.2 (20CRV2C) to +0.3 • C decade −1 (ERA20C) at 700 hPa, and reaching +0.3 • C decade −1 (ERA20C and 20CRV2C) at 500 hPa. However, ERA20C and 20CRV2C also show an increase in HUS, while ERA40/ERA-Interim and NCEP/NCAR-v1 show a decrease. This increase in HUS is significant in ERA20C at both pressure levels, and is valued at +0.054 g kg −1 decade −1 at 700 hPa and +0.022 g kg −1 decade −1 at 500 hPa. In addition, compared with ERA40/ERA-Interim and NCEP/NCAR-v1, ERA20C is dryer and colder at both investigated pressure levels, whereas 20CRV2C is warmer and wetter ( Figure S8, Supplementary Materials). These trends suggest that decreasing LCC and MCC during the brightening period could result from increasing TA and decreasing HUS generating a decrease in HUR, which could have therefore inhibited cloud development. The initial changes in TA and HUS could be explained by changes in the atmospheric circulation. Eastman and Warren [62] attribute the declining cloud cover observed over Western Europe to a poleward shift of the polar jet stream since the 1970s. Bender et al. [66] highlighted a poleward shift of extratropical storm tracks (which induces a poleward shift of the polar jet stream) for 1983-2008 as a consequence of global warming. This poleward shift of extratropical storm tracks is accompanied by a reduction in total cloud cover, leading to a positive cloud feedback [66]. The conclusions of these studies are not in opposition with our hypothesis, since increasing TA during the brightening period (1980-2010) could result from more frequent tropical air advections over Western Europe. In addition, we know that the brightening period was characterised by a strong decrease in the aerosol loading in the troposphere [59,60]. This decline in aerosol loading could be responsible for the decrease in LCC and MCC during the brightening period by inhibiting cloud development. Then, the decrease in cloud cover could have generated more E g↓ , and therefore increasing TA, which could generate an decrease in HUR.
In summary, there are two hypotheses to explain the changes in cloud cover during the dimming and the brightening periods: the "atmospheric circulation hypothesis" and the "aerosol cloud interactions hypothesis". On the one hand, atmospheric circulations changes could have modified TA and HUS, and therefore HUR, with consequences on cloud development and later E g↓ (atmospheric circulation hypothesis). On the other hand, aerosol-cloud interactions could have modified cloud development, and therefore E g↓ , with consequences on TA and HUR. We cannot settle on a single hypothesis to explain LCC and MCC changes (especially since statistics are rarely significant), but they are most probably due to the combination of several factors. There are indeed proofs for both hypotheses, especially during the brightening period [59,62,66].
It should be noted that the similar trend values found in ERA40/ERA-Interim and NCEP/NCAR-v1 at 700 hPa and 500 hPa confirm that the apparent artefact found in the TCC computed by MAR-NCEP1 was the result (Section 3.1) of biases in the forcing fields of MAR in the high troposphere (above 440 hPa).

Conclusions
The main objective of this paper is to assess the ability of the MAR RCM to reproduce observed changes in E g↓ in Belgium with respect to reanalyses and observations. To this end, MAR was forced by different reanalyses: ERA-Interim (1979-2010) completed by ERA40 (1958ERA40 ( -1978, NCEP/NCAR-v1 (1958-2010), ERA20C (1958ERA20C ( -2010, and 20CRV2C . This study focusses on the 1959-2010 period, which is the longest common period covered by all of the reanalyses used as MAR forcing. In addition, we consider two distinct periods in our analysis: 1959-1979 (dimming) and 1980-2010 (brightening). Measurements of E g↓ from GEBA and RMI, as well as cloud cover observations from Belgocontrol and RMIB, were used for the evaluation of MAR and the forcing reanalyses.
The comparison of E g↓ and TCC modelled by MAR and reanalyses with measurements performed at several ground stations show that MAR allows us to largely reduce the mean biases present in the reanalyses. However, MAR does not significantly improve R and CRMSE, as the interrannual variability of RCMs is prescribed by their forcings. The comparison also shows that the biases are the largest in summer for both E g↓ and TCC, and that MAR generally overestimates E g↓ , while it underestimates TCC. This suggests that the biases in E g↓ result in large part from inaccuracies in TCC. Improvements in cloud physics are still required, despite the modification performed in the new version of MAR (MAR v3.8). The comparison with forcing reanalyses finally suggests that the cloud optical depth is larger in MAR than in the reanalyses, as the reanalyses show larger biases in E g↓ , while MAR and reanalyses exhibit similar biases in TCC.
The analysis of the linear trends in annual and seasonal E g↓ and low, medium, and high clouds brings to light that only MAR-ERA successfully simulates consistent trends in E g↓ , which is probably because ERA40/ERA-interim has a better spatial resolution and assimilates more observations than the other reanalyses used in this study. Neither MAR-ERA-20C nor MAR-20CRV2C shows any of these trends, probably because both ERA-20C and 20CRV2C reanalyses assimilate surface observations only. MAR-NCEP1, for its part, overestimates the trend computed for the brightening period because of biases in the forcing fields in the high troposphere, especially before the 1990s. Our results show that annual E g↓ trends are mainly driven by E g↓ changes in spring and summer. They also suggest that in Belgium, the aerosol-radiation interactions cannot be neglected, especially for the dimming period, because trends in observed TCC are non-significant, while trends in observed E g↓ are significant. The increase in E g↓ that was observed in Belgium since the 1980s and especially since the 2000s could mainly be explained by a decrease in the low and medium cloud cover strengthening the effect of the decrease in aerosol loading on E g↓ that has been observed in Europe since the 1980s. The origin of these changes in cloudiness is not clear, and could result from changes in both atmospheric circulation and aerosol-cloud interactions. Further analyses are therefore required. First of all, as this trend is not simulated when MAR is forced by 20CRV2C and ERA-20C, comparing the MAR forcing fields (temperature, humidity, and wind) from these 20th century reanalyses with ERA-Interim will help detect the meteorological variables that are responsible for this decrease in cloudiness. Sensitivity experiments with MAR could be performed using variables from two different reanalyses, e.g., wind from ERA-Interim, and other fields from ERA-20C. Afterwards, a circulation type classification (as the one developed by Belleflamme et al. [67]) could be used to associate the daily cloudiness to some types of general circulation. The temporal evolution of these circulation types could then be studied as was done in Bednorz et al. [68] and Lupikaszca and Lipinski [69]. Finally, it might also be interesting to verify whether such a decrease in cloudiness would be projected if MAR was forced by GCM-based future projections. This would suggest that such a change could be a consequence of global warming.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4433/9/7/262/s1. Figure S1. Modelled E g↓ vs. observed E g↓ at Melle for 1967-2010: (a) annual and seasonal correlation (R), (b) mean bias (MB), (c) centred root mean squared error (CMRSE) and (d) standard deviation (STD). It should be noted that ERA40 covers 1967-1978 while ERA-interim covers 1979-2010; Table S1. Seasonal trends and their uncertainty range of Snedecor for the 95% confidence interval computed computed from measurements performed at Saint-Hubert. The values in bold indicate that the trends are statistically significant at the 95% confidence level, namely that they are larger than the associated uncertainty range. The values between brackets indicate that they are computed for 1966-1979; Figure S2.