An integrated and homogenized global surface solar radiation dataset and its reconstruction based on a convolutional neural network approach

Surface solar radiation (SSR) is an essential factor in the flow of surface energy, enabling accurate capturing of long-term climate change and understanding of the energy balance of Earth's atmosphere system. However, the long-term trend estimation of SSR is subject to significant uncertainties due to the temporal inhomogeneity and the uneven spatial distribution of in situ observations. This paper develops an observational integrated and homogenized global terrestrial (except for Antarctica) station SSR dataset (SSRIH_station) by integrating all available SSR observations, including the existing homogenized SSR results. The series is then interpolated in order to obtain a 5^∘ × 5^∘ resolution gridded dataset (SSRIH_grid). On this basis, we further reconstruct a long-term (1955–2018) global land (except for Antarctica) SSR anomaly dataset with a 5^∘ × 2.5^∘ resolution (SSRIH_20CR) by training improved partial convolutional neural network deep-learning methods based on 20th Century Reanalysis version 3 (20CRv3). Based on this, we analysed the global land- (except for Antarctica) and regional-scale SSR trends and spatiotemporal variations. The reconstruction results reflect the distribution of SSR anomalies and have high reliability in filling and reconstructing the missing values. At the global land (except for Antarctica) scale, the decreasing trend of the SSRIH_20CR (−1.276 ± 0.205 W m⁻² per decade) is smaller than the trend of the SSRIH_grid (−1.776 ± 0.230 W m⁻² per decade) from 1955 to 1991. The trend of the SSRIH_20CR (0.697 ± 0.359 W m⁻² per decade) from 1991 to 2018 is also marginally lower than that of the SSRIH_grid (0.851 ± 0.410 W m⁻² per decade). At the regional scale, the difference between the SSRIH_20CR and SSRIH_grid is more significant in years and areas with insufficient coverage. Asia, Africa, Europe and North America cause the global dimming of the SSRIH_20CR, while Europe and North America drive the global brightening of the SSRIH_20CR. Spatial sampling inadequacies have largely contributed to a bias in the long-term variation of global and regional SSR. This paper's homogenized gridded dataset and the Artificial Intelligence reconstruction gridded dataset (Jiao and Li, 2023) are both available at https://doi.org/10.6084/m9.figshare.21625079.v1.

1 Introduction

Energy flows at the Earth's surface play an essential role in climate change and human activity and link to physical processes such as global warming, glacier retreating, the hydrological cycle and the carbon budget (Hoskins and Valdes, 1990; Peixoto et al., 1992; Trenberth and Fasullo, 2013; Wild, 2012). As a critical factor characterizing surface energy flows, surface solar radiation (SSR) largely determines the climatic conditions and ecological environment in which we live. Therefore, a more accurate and comprehensive analysis of the SSR fluxes will help better understand the Earth's atmospheric system. In situ observations provide the most accurate baseline data for measuring SSR. They allow for the first time the detection of decadal changes in SSR known as “dimming and brightening” (Wild et al., 2005), especially considering that they cover a longer period concerning another type of data, e.g. satellite data (Pfeifroth et al., 2018). Even observational data often have uneven distribution and missing data with respect to the satellite data, especially in areas with complex orography (Manara et al., 2020).

The sources of in situ SSR observations are mainly collected from the Global Energy Balance Archive (GEBA) (Wild et al., 2017) and the World Radiation Data Centre (WRDC) (Tsvetkov et al., 1995). Furthermore, other SSR station series are obtained from the high-quality Baseline Surface Radiation Network (BSRN) (Driemel et al., 2018) and the data centres of individual national hydrometeorological services. However, two issues still need to be addressed: (1) the inhomogeneity of station data resulting from station relocations and instrumentation changes severely impacting the climate change assessment. For the regions with a relatively high density of stations, like Europe (Manara et al., 2019, 2016; Sanchez-Lorenzo et al., 2013a, b, 2015), Japan (Ma et al., 2022) and China (Ju et al., 2006; Wang, 2014; Wang et al., 2015; Wang and Wild, 2016; S. Yang et al., 2018; You et al., 2013), much previous work has redefined the degree and timing of dimming and brightening by addressing the inhomogeneity of the SSR data series. For example, in Spain, the average annual homogenized SSR series has a significant increasing trend (+3.9 W m⁻² per decade) during the 1985–2010 period (Sanchez-Lorenzo et al., 2013a). The period of dimming observed in Italy's homogenized SSR series is not apparent in the 1960s and early 1970s, when the raw series (non-homogenized) are taken into account (Manara et al., 2016). The direct measurements of SSR show a level trend from 1961 to 2014 over Japan, while their homogenization series display a decreasing trend (0.8–1.6 W m⁻² per decade) (Ma et al., 2022). In China, homogenization largely eliminated the dramatic non-climatic rise of the early 1990s and also reduced the increasing trend from 1990 to 2016 (S. Yang et al., 2018). However, most of the research was still limited to regional scales. (2) There is the issue of limited spatial sampling of long observational stations and their uneven distribution, especially over areas with complex orography. Considerable efforts have been devoted to filling in or interpolating the missing values in climate datasets (“spatial analysis”) (Collins, 1996; Erxleben et al., 2002; Scudiero et al., 2016). The traditional spatial interpolation methods commonly used include inverse distance weighting (Fisher et al., 1993; Shepard, 1968), kriging (Krige, 1951) and thin-plate splines (Bookstein, 1989). Since the 1980s, physical parametric interpolation (Feng and Wang, 2021; Tang et al., 2019) and Bayesian fusion schemes (Aguiar et al., 2015) based on multi-source observational data have been widely used with the emergence of highly accurate and relatively precise satellite data. However, the resulting fusion datasets cover too short a period to investigate their decadal and multi-decadal variations and to study the underlying causes. The spatial, temporal and spectral coverage of a single satellite is limited, and multiple satellite data are therefore often used in tandem with each other; however, such a discontinuity in time and space can introduce inhomogeneity into a dataset (Evan et al., 2007; Feng and Wang, 2021; Shao et al., 2022). Reanalysis products are an important complement containing long-term SSR data and therefore have been widely used in climate studies (Huang et al., 2018; Jiao et al., 2022; Urraca et al., 2018; C. Zhou et al., 2018, 2017) due to the dynamically consistent and spatiotemporally complete atmospheric fields with high resolution and open access to data. However, existing studies have shown that reanalysis products generally overestimate multi-year mean SSR values compared to observations over land (He et al., 2021). With the continuous development of climate system simulations, model data from the Coupled Model International Project (CMIP) have become an important resource for conducting climate change research (Gates et al., 1999; Zhou et al., 2019). Previous studies have shown that the models used in CMIP6 overestimate the global mean SSR (He et al., 2023; Jiao et al., 2022; Wild, 2020). The rise of deep-learning and big-data techniques has brought about an explosion of artificial intelligence (AI). Machine learning is increasingly being used in spatial interpolation, such as the spatial reconstruction of surface temperature datasets (Huang et al., 2022; Kadow et al., 2020) or the spatial and temporal reconstruction of turbulence resolution (Fukami et al., 2021). Furthermore, it shows high accuracy and low uncertainty in reproducing and predicting SSR (Leirvik and Yuan, 2021; Tang et al., 2016; L. Yang et al., 2018; Yuan et al., 2021). However, long-term homogenized SSR datasets with global terrestrial coverage have yet to be developed, resulting in significant uncertainties in assessing global SSR variation (Jiao et al., 2022).

Therefore, developing a more homogeneous and comprehensive global long-term SSR climatic dataset that provides a better benchmark for observational constraints on the global surface energy balance and budget remains a valuable and challenging task. This paper first homogenizes and grids the most extensive collection of available global SSR station observations. Then, the missing grid boxes and years are spatially interpolated using a convolutional neural network (CNN) approach to obtain a globally covered land surface SSR anomaly dataset. Finally, the reconstructed datasets are initially analysed and evaluated. Thus, the paper is divided into seven main sections. The data resources are introduced in Sect. 2. Section 3 presents the data homogenization and the CNN model reconstruction methods. The data homogenization and verification are shown in Sect. 4. Section 5 gives the AI reconstruction results. Section 6 is the availability of the datasets. Conclusions are provided at the end of the paper.

2 Data

Nine SSR datasets are collected to derive the global SSR variable. In particular, six datasets contain data from observational stations (Sect. 2.1): two global ground-based measurement datasets (GEBA, WRDC) and four homogenized products at the regional and country levels (Europe, China, Japan and Italy). Three of the adopted datasets are reanalysis data (Sect. 2.2.1): fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis (ERA5), 20th Century Reanalysis version 3 (20CRv3) data and the Coupled Model Intercomparison Project Phase 6 (CMIP6) historical simulation output (125). Specifically, the ERA5 data are used to fill the data over oceans and Antarctica (Sect. 3.2.1), and 20CRv3 data and CMIP6 simulations are used for AI model training (Sect. 5.1) and reconstruction. All are listed in Table 1.

Table 1 List of information on the various types of data used in this paper.

2.1 In situ observational data

2.1.1 Global datasets

There are two main sources of raw SSR data (see Table 1): the ETH Zurich GEBA with monthly data from 2445 globally distributed stations, starting from 1922 until 2020, and the WRDC dataset with monthly globally distributed data from 1136 stations since 1964. The first one is available for download at https://geba.ethz.ch (last access: 2 July 2022). The second one published the first SSR radiation balance data in 1965, and its publication has been issued four times a year since 1993 and is available for download at http://wrdc.mgo.rssi.ru/ (last access: July 2021).

2.1.2 National (regional) homogenized station datasets

(1) Chinese homogenized SSR dataset

The China Meteorological Radiation Fundamental Elements Monthly Value Data Set was downloaded from http://www.nmic.cn (last access: September 2022). The homogenized SSR dataset in China is released by the National Meteorological Information Centre (NMIC) of the China Meteorological Administration (CMA) (S. Yang et a., 2018). The data are available for the period between January 1950 and December 2014, and the follow-up data are extended with raw observations from the NMIC. They used the sunshine duration (SSD) data from nearby stations to construct an arguably better reference to identify inhomogeneities in the SSR data. Then, a combined metadata and maximum-penalty t-test (PMT) method was used to detect the change points. Finally, they were adjusted by a quantile-matching (QM) algorithm (Wang and Feng, 2013). The final homogenized SSR station dataset was converted to gridded data using the first difference method (FDM, Peterson et al., 1998) and is available for download at http://www.nmic.cn (last access: September 2022).

(2) Japanese homogenized SSR dataset

Ma et al. (2022) released a Japanese SSR homogenized dataset in 2022 spanning the period between 1870 and 2015. First, they homogenized SSD based on PMF (penalized maximal F test) and QM algorithms. They then used the homogenized SSD from the previous step as a reference series, combined with metadata and PMT, to detect change points. Finally, they adjusted the change points by the QM algorithm. For more details on data descriptions, the adopted methodology and downloading data, see https://data.tpdc.ac.cn/en/data/45d73756-3f5a-4d27-82a4-952e268c20e8/ (last access: March 2022).

(3) European homogenized SSR data

A homogenized dataset of European SSR stations was developed by Sanchez-Lorenzo et al. (2015) and is currently available for full public download at https://doi.org/10.1002/2015JD023321. They selected the 56 longest central European SSR series available in the GEBA dataset with data for the period between 1922 and 2012. They adjusted them to ensure temporal homogeneity, homogenizing the data with the standard normal homogeneity test (Alexandersson, 1986) and the Craddock test (Craddock, 1979).

(4) Italian homogenized SSR dataset

The Italian homogenized SSR datasets are those published by Manara et al. (2019, 2016). As candidate stations to use as reference series, they selected the 10 series located in the same area of the series to be tested, and that series correlates well with the test one. In particular, they tested the change points with the Craddock test (Manara et al., 2017), and when a break is identified by more than one reference series, the preceding portion of the series is corrected, leaving the most recent portion unchanged. In this way, the SSR stations were homogenized, and then the missing values were interpolated.

2.2 Other datasets

2.2.1 Reanalysis

ERA5 can be used to fill in SSR data from the oceans and Antarctica and carry out the global reconstruction, taking into account its high spatial resolution and the reliable performance of SSR (Jiao et al., 2022; Liang et al., 2022). After the reconstruction, we removed the data for the ocean reanalysis and maintained the data only in the land area (except for Antarctica). In addition, two SSR data products (20CRv3, CMIP6) are used to train AI models. These are the following.

1. ERA5 (space-filling data). ERA5 is the fifth generation of the European Centre for Medium-Range Forecasts reanalysis product, which currently publishes data from 1950 to the present (Hersbach et al., 2020). In addition, ERA5 has an hourly output and an uncertainty estimate from the ensemble. The data are based on the Integrated Forecasting Model Cy41r2 run in 2016, which contains a 4D-Var assimilation scheme. In ERA5, SSR is obtained from a rapid radiation transfer model (RRTM) (Mlawer et al., 1997). The present study utilizes monthly SSR data for the period 1955–2018 from ERA5 with a resolution of 0.25^∘ × 0.25^∘. They can be downloaded at https://cds.climate.copernicus.eu (last access: July 2022).

2. 20CRv3 (data for AI model training). The 20CR project is an effort led by the NOAA's Physical Sciences Laboratory and CIRES at the University of Colorado, supported by the Department of Energy, to produce reanalysis datasets spanning the entire 20th century and much of the 19th century (Slivinski et al., 2019). 20CR provides a comprehensive global atmospheric circulation dataset from 1850 to 2015. Its chief motivation is to provide an observational validation dataset, with quantified uncertainties, for assessing climate model simulations of the 20th century. 20CR uses an ensemble filter data assimilation method which directly estimates the most likely state of the global atmosphere every 3 h and estimates the uncertainty in that analysis. The most recent version of this reanalysis, 20CRv3, provides 8-times daily estimates of global tropospheric variability across 75 km grids, spanning 1836 to 2015 (with an experimental extension from 1806 to 1835). The present study uses monthly SSR data of 20CRv3 (NOAA/CIRES/DOE 20CR, 80 members) from 1955 to 2015. We selected all 80 members of the 20CR as input (1 for evaluation and to test reconstruction, the other 79 for training the CNN model). The SSR of 20CRv3 has a spatial resolution of 0.7^∘ × 0.7^∘. The download is available at https://portal.nersc.gov/archive/home/projects/incite11/ (last access: May 2022).

2.2.2 CMIP6 model output

• 3.

  CMIP6 model output (data for AI model training). The Coupled Model Intercomparison Project, driven by the World Climate Research Program, is now in its sixth phase. Specifically, CMIP6 is considered the current state-of-the-art way of producing future climate simulations, including predicting future SSR based on different climate scenarios (W. Zhou et al., 2018). It provides an important resource for studying current and future climate change (Eyring et al., 2016). The historical simulations of CMIP6 are designed to reproduce observed climate and climate change constrained by radiative forcing. CMIP6 historical simulation spans between 1850 and 2014. In this study, we selected 125 members out of a total of 507 members from several CMIP6 large-ensemble models (with more than 10 realizations and runs) with high correlation coefficients with observations as input to train and validate the CNN model (1 for evaluation and to test reconstruction, the other 124 for training the CNN model). We selected the monthly downward shortwave radiation from 1955 to 2014 (see Table S1 in the Supplement). The data can be downloaded at https://esgf-node.llnl.gov/search/cmip6 (last access: July 2022).

3 Methods

3.1 Data quality control (QC) and homogenization

The SSR data homogenization method is only applied to the two non-homogenized in situ observation datasets (GEBA and WRDC). The QC and homogenization flowchart (Fig. 1) is divided into three steps: (1) QC; (2) homogenization; and (3) integration and consolidation.

Figure 1 Flowchart of quality control (QC) (first step), homogenization (second step) and integration (third step).

3.1.1 QC

The QC of SSR data includes the following steps.

1. Simple integration is integration of the GEBA (2445) and WRDC (1136) datasets, removing stations with no data and leaving 2681 stations.

2. Removing duplicate stations. (a) For stations with similar latitude and longitude, we consider two stations with totally identical latitude and longitude to be the same station. (b) For stations less than 10 km apart, we averaged the duplicate stations in these a and b cases. (c) For special duplicate stations, we stitched together data of the duplicate stations based on metadata from the CMA.

3. Remove stations, years or months for which a climatic analysis cannot be established: we remove stations with records of less than 10 years and values more than 3 times (3σ criterion, Olanow and Koller, 1998) the standard deviation of the SSR anomalies.

4. Candidate stations (487) with a record length greater than 15 years in the period 1971–2000 are selected. We added stations (715) with more than 10 years of SSR records to increase the number of available stations for a better homogenization of the candidate stations (Fig. 2).

Figure 2 Spatial distribution of candidate stations (“*”) and added stations (“+”). The different colour bars represent the length of the station record in months.

3.1.2 Station series homogenization

This paper uses the RHtestV4 software package to test and adjust the SSR station data for homogeneity (http://etccdi.pacificclimate.org/software.shtmlm last access: July 2021) (Wang and Feng, 2013). The package is based on the empirical penalty functions PMF (Wang, 2008a) and PMT (Wang, 2008b; Wang et al., 2007) for the homogenization test. It takes into account the lag-1 autocorrelation of the time series. It embeds a multiple linear regression algorithm to significantly reduce the problem of an unbalanced distribution of pseudo-identification rates and test efficacy. Also, RHtestV4 uses the QM algorithm (Vincent et al., 2012; Wang et al., 2010) and mean adjustments to adjust the identified change points.

The specific steps are as follows.

• 1.

  Building the reference series

  • a.

    We processed the data from all station series (715) into the annual first difference (FD) series e_i (Eq. 1) (Peterson et al., 1998).

  • b.

    We calculated the correlation of the annual FD series between the series from the potential reference pool and the candidate stations.

  • c.

    We calculated the distance between the potential reference pool stations and candidate stations.

  • d.

    We selected potential stations according to the correlation coefficient (CC > = 0.6) between the series from potential reference pool and candidate stations. The potential stations also satisfy the limits in distances (< = 500 km) between the potential pool stations and candidate stations.

  • e.

    We obtain the reference FD series (Re) based on the m potential reference series (Pe_i) and the CCs (c_i) between the potential reference series (Pe_i) and candidate station series (Eq. 2).

  • f.

    The synthesized reference FD series (Re) (Eq. 2), plus the average of all potential reference series ($\stackrel{\mathrm{‾}}{R}$), yields the final annual reference series (R) (Eq. 3).

    $$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{Pe}}_i=x_i-x_{i+\mathrm{1}}i=\mathrm{1},\mathrm{2},\mathrm{\dots },n-\mathrm{1}\text{(2)}& {\displaystyle }\mathrm{Re}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m{\mathrm{Pe}}_i\cdot c_i^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^m{c}_i^{\mathrm{2}}}}\text{(3)}& {\displaystyle }R=\mathrm{Re}+\stackrel{\mathrm{‾}}{R}\end{array}$$

    x_i is the raw observational station SSR in the year i, Re is the final reference series, Pe_i is the potential reference series, and c_i is the CC between the potential reference series and the candidate station series.

• 2.

  Testing and adjusting the candidate series

  The homogenization test algorithm used in this paper is the PMT. This method is a reference series-dependent test for a normalized candidate series. It assumes that the linear trend of the time series is zero and uses the degree of mean deviation at different points in the series to find change points. Furthermore, it eliminates the effect of different sample lengths on the test results. At the same time, the method introduces an empirical penalty factor, which effectively improves detection. We used the PMT to test the homogeneity of the candidate series based on the reference series established in (1). We then adjusted the statistically significant (p>0.05) change points obtained using the mean adjustment method (p>0.05). We homogenized the monthly series for 66 stations (see Fig. S1 in the Supplement).

3.1.3 Integration and consolidation

As can be seen from Fig. 1, the candidate stations (487) are relatively sparse. To better adapt deep-learning methods for the dataset reconstruction later, we adjusted, added and integrated station series based on the results of homogenized data from other scholars. (1) We added stations with more than 10-year overall (1955–2018) records but no more than 15 years during the 1971–2000 period and removed those stations that were clearly inhomogeneous (25) and some years of station (3). (2) We subsequently integrated monthly SSR series for 116 stations based on the results of homogenization from other scholars: China (56), Japan (8), Europe (2) and Italy (50). After the above steps, we ended up with a homogenized dataset containing 944 stations (Fig. 3). The details of the processing and classification are shown in Table S2 (see the Supplement).

Figure 3 Spatial distribution of stations after homogenization (unit months). Different colours represent the length of the station records in months.

3.2 CNN model reconstruction methods

The CNN deep-learning model network architecture uses a U-shaped structure similar to U-Net (Ronneberger et al., 2015). The advantage of using this model is that (1) both high- and low-frequency information of the picture can be retained, and when reconstructing the SSR data, not only will the grid point information close to the missing measurement point be considered, but information from more distant locations will be too (which may be remotely correlated with that missing measurement point). (2) This makes the model convergence faster and more economical in terms of computational resources. The upper part of the U-shaped structure, which has no downsamples or a low number of downsamples, represents the high-frequency information of the graph. These sections contain much of the detail in the graph, and the relationships between similar grid points are conveyed by this section. The lower half of the U-shaped structure is downsampled more often and represents the lower-frequency information of the graph. The global radiation of a wide range of undulations is transmitted by it, and then the information at the various levels of the U-shaped structure is connected and transmitted through the skip connection, allowing the whole network to remember all the information of the picture very well. The model uses nearest-neighbour upsampling in the decoding phase, and the skip links will concatenate two feature maps and two masks as the feature and mask inputs for the next part of the convolution layer. The input to the last part of the convolution layer will contain the original input image concatenated with the holes and the original mask, allowing the model to replicate the gap-free pixels. The complex and variable nature of the sea–land boundary then has a significant impact on the reconstruction when we reconstruct the global land SSR data. Therefore, we use partial convolution at the image boundaries with a suitable image padding, ensuring that the padding content at the image boundaries is not affected by values outside the image. The deep-learning models' convolutional layers and loss functions are described in the Supplement.

We further reconstruct a long-term (1955–2018) global SSR anomaly dataset (SSRIH_20CR) by using improved partial CNN deep-learning methods based on a “perfect” dataset. A CNN consists of three parts: a convolutional layer to reduce the number of weights by extracting local features, a pooling layer to reduce peacekeeping and prevent overfitting, and a fully connected layer to output the desired result. In this paper, a modified CNN is used to model the reconstruction of the SSR data, with the convolutional layer replaced by a partial convolution method and mask update. This method is the latest in image restoration effects and can restore irregular holes, an advantage over other image restoration methods that can only restore rectangular holes. Therefore, this paper uses the modified CNN model (Kadow et al., 2020) to recover the missing part of the global terrestrial SSR (except Antarctica). The specific reconstruction steps and processes are as in Fig. 4.

Figure 4 Flowchart of AI reconstruction.

3.2.1 Data pre-processing

The homogenized station data are converted to grid box anomalies using the climate anomalies method (CAM) (Jones et al., 2001). CAM is a commonly used method for converting station anomaly data to gridded data. We divide all global areas into a 5^∘ × 5^∘ grid, after which we calculate the SSR anomalies (relative to 1923–2020) within the grid box by averaging the anomalies of all the stations (at least one station in it). If more than one site exists in the same grid box, the record length of this grid box is the total length of all sites in that grid box. Finally, we removed the values that were more than 3 times the standard deviation of the SSR anomaly time series after gridding. SSRs are all processed as daily average anomalies, i.e. monthly anomalies divided by 30 (each month is approximated as 30 d). We multiply all the values by 30 again when the reconstruction is complete. The global land (except for Antarctica) distribution and coverage of SSRs after gridding are shown in Fig. 5a, b.

Figure 5 (a) Spatial distribution of 5^∘ × 5^∘ grid boxes (SSRIH_grid) obtained by interpolating the homogenized global land (except for Antarctica) SSR series. The different colours represent the length (the sum of all the records) of the station record in unit years. (b) Grid box coverage for the homogenized global land (except for Antarctica) SSR (SSRIH_grid) except for Antarctica.

As seen in Fig. 5a, the SSR is spatially sparsely distributed across South America and Africa. As shown in Fig. 5b, SSR coverage increased yearly from 1950 until the mid-1970s, when it slowly decreased. In 2013, the coverage rate decreased sharply due to untimely data submission. Considering the SSR coverage above, we only kept the years (1955–2018) with data coverage of more than 8 % of global land (except for Antarctica) areas.

Comparisons show that ERA5 has a high spatial resolution and relatively reliable performance in the temporal variations and long-term trends (Liang et al., 2022; Jiao et al., 2022). To obtain a higher data coverage and ensure that the AI model runs well, we used the ERA5 to fill the SSR of the homogenized global gridded SSR in the Antarctic and ocean areas. However, if we use the SSR of ERA5 to directly fill the SSR of the homogenized global gridded SSR (SSRIH_grid) in the Antarctic and in the ocean areas, then the relatively weaker ocean SSR variations (variabilities, decadal changes, trends) from ERA5 will inevitably introduce certain systematic biases in land SSR reconstruction due to the SSRs having the lower coverage on the land. Therefore, we designed an algorithm to avoid excessive diffusion of SSR system bias in terrestrial areas: we first calculated the ratios ${\mathit{\gamma }}_i(i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },n)$ between the SSR from ERA5 and from SSRIH_grid on the land in all n years. For a single grid box, the γ_i have small changes and are regarded as a constant γ_median (Eq. 4), and the γ_median vary by latitude and longitude in both the marine and land areas. We then extrapolated the γ_median for all the grid boxes along the land and sea boundaries. If there is no observation there, then the adjacent ocean ERA5 SSR is used to take its place after it is adjusted according to the differences between the SSR variations (represented by the linear trends) for the different underlying surfaces (Eq. 5):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{median}}=\mathrm{median}\left({\displaystyle \frac{{\mathrm{OBS}}_{i\mathit{\_}\mathrm{land}}}{\mathrm{ERA}{\mathrm{5}}_{i\mathit{\_}\mathrm{land}}}}\right),\text{(5)}& {\displaystyle }\begin{array}{rl}{\mathrm{OBS}}_{i\mathit{\_}\mathrm{O}\mathit{\&}\mathrm{L}}\left(\mathrm{land}\right)& =\mathrm{ERA}{\mathrm{5}}_{i\mathit{\_}\mathrm{O}\mathit{\&}\mathrm{L}}\left(\mathrm{ocean}\right)\cdot {\mathit{\gamma }}_{\mathrm{median}}\cdot {\displaystyle \frac{T_{\mathrm{O}}}{T_{\mathrm{L}}}}\\ & i=\mathrm{1},\mathrm{2},\mathrm{3},n.\end{array}\end{array}$$

γ_median is the median value of the ratios of observation (OBS) and ERA5 land SSR series. OBS_{i_land} is the land SSR for the year i from the SSRIH_grid in a single grid. ERA5_{i_land} is the land SSR for the year i from ERA5 in a single grid. OBS_{i_O&L}(land) is the land SSR along the sea–land boundary (land) for the year i from the SSRIH_grid. ERA5_{i_O&L}(ocean) is the ocean SSR along the sea–land boundary for the year i from ERA5. T_O is the trend of the ERA5 SSR in ocean areas in all n years, and T_L is the trend of the ERA5 SSR in areas in all n years.

3.2.2 AI model reconstruction

We use a server (configured with processor Intel (R) Core (TM) i7-8700 CPU @ 3.20 GHz 3.19 GHz, RAM 32G, 64-bit OS, GPU model 516.94, NVIDIA GeForce 1080T version, Python 3.9.12 64-bit, CUDA 10.1) for AI model training. The specific training steps are as follows.

1. A total of 768 missing-value masks (monthly masks between 1955 and 2018) were prepared for training and validation using “1” for existing and “0” for missing values.

2. The 20CRv3–CMIP6 training set (monthly values between 1955 and 2015/2014) and missing-value masks are fed into the 20CR-AI and CMIP6-AI model for training.

3. We perform 1 500 000 training sessions with an interval of 10 000 sessions for the training output model.

Afterwards, the two AI models are validated against the root mean squared error (RMSE) and CCs of the reconstructed SSRs (SSR${}_{\mathrm{20}\mathrm{CR}}/$SSR_CMIP6). The validation set SSRs, and the optimal number of training cycles is 1 100 000 (see Figs. S2, S3 and S4 in the Supplement). The initial hyperparameters of the model are set as follows: a learning rate of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ and learning fine-tuning of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$. First, we set the batch size to 16 in the first 500 000 iterations and fine-tune it to 18 in the last 10 000 000 iterations, for a total of 1 500 000 iterations, to suppress the overfitting phenomenon generated during the training process. We validate the model every 10 000 times and with early stopping if the validation shows a decreasing trend, and the final number of training times used is 1 100 000. Second, L2 (ridge regression) regularization is also added to regulate the loss function (see Eq. S9 in the Supplement).

The training result models generated by the different AI models are obtained separately for the different training sets. The model is first used to reconstruct a reanalysis validation set with the same missing-value mask as the original observation dataset. This is followed by a validation of the reconstruction against the original reanalysis dataset (calculation of CC and RMSE) to understand the discrepancies in the model reconstruction.

4 Data homogenization and verification

We homogenized the original monthly station or gridded SSR time series (SSRIH_station or SSRIH_grid) using the method in Sect. 3.1.2. We selected six continental regions, excluding Antarctica and the Arctic, from the eight continents of the world defined by Xu et al. (2018) (Asia, Africa, South America, Europe, North America, Australia, Antarctica and the Arctic). The decreasing trend of the SSRIH_grid is consistent with the original gridded SSR series (SSRI_grid) during 1955–1991, while the increasing trend during 1991–2018 is weaker. At the regional scale, the SSRIH_grid has a generally similar variation to the SSRI_grid, and the SSRIH_grid is usually more representative of climate change than the SSRI_grid at individual stations.

Figure S5 (see the Supplement) illustrates the long-term variations of global (Fig. S5a in the Supplement) and continental land SSR (Fig. S5b in the Supplement) from the SSRI_grid and SSRIH_grid (except for Antarctica) during 1955–2018. The most prominent change revolves around the adjustment around 1992: the SSR anomalies were systematically adjusted upward from 1987 to 1992, while the SSR anomalies were systematically adjusted downward from 1993 onwards. Thus, there is a significant decreasing trend for both the global land SSRI_grid (−1.995 ± 0.251 W m⁻² per decade) and global land SSRIH_grid (−1.776 ± 0.230 W m⁻² per decade) (except for Antarctica) from 1955 to 1991, while the increasing trend of the global land SSRIH_grid from 1991 to 2018 is 0.851 ± 0.410 W m⁻² per decade, slightly smaller than the increasing trend of the SSRI_grid (0.999 ± 0.504 W m⁻² per decade). It is worth noting that 1992 happened to be the second year of the eruption of Mount Pinatubo, and the homogenized SSR data integrated in this paper may be affected by this event. However, overall, the homogenization also has limited effects on the global SSR variations from Fig. S5 (see the Supplement), which is consistent with the influence of data homogenization on a wide range of surface air temperatures (Brohan et al., 2006; Xu et al., 2013).

On the regional scale, the differences between the SSRIH_grid and SSRI_grid are more pronounced in Asia and Europe (see Fig. S5b in the Supplement). Asia's homogenized SSR shows that the regional average SSR has been declining significantly over the period 1958–1990; this dimming trend mostly diminished over the period 1991–2005 and was replaced by a brightening trend in the recent decade. The SSRIH_grid in Asia is higher than the SSRI_grid from 1985 to 1990 and lower than the SSRI_grid from 2012 to 2015. The SSRIH_grid shows a more moderate short-term increase in Europe from 1960 to 1980. Note also that the Australian raw data prior to 1988 were artificially detrended because at the time the Australia Weather Service was afraid that the instruments would drift. Therefore, they detrended them and unfortunately did not store the raw data, and the SSR evolution in Australia is artificial with no trend (Wild et al., 2005). In addition, the SSRI_station and SSRIH_station comparisons for all 66 stations are shown in Fig. S1 (see the Supplement).

5 AI reconstruction and comparison

5.1 Training of the AI model

We produce two (20CRv3 and CMIP6) separate training and validation sets: we select the first member data of the reanalysis data and the model data, respectively, as the validation set, and the remaining 79 (124) ensemble members as the training sets, where each ensemble member included 732 (720) months of SSR data. Each validation set included 732 (720) samples, while the training sets contained 57 828 (89 280) ensemble members. All the above data, including the in situ observations, are then resampled to monthly anomalies of 5^∘ × 2.5^∘. We reconstruct the SSR of 20CRv3 and CMIP6 with missing values based on the 20CRv3 and CMIP6 datasets using the method in Sect. 3.2 and obtain two reconstructions, SSR_20CR and SSR_CMIP6, respectively. The SSR of 20CRv3 and CMIP6 with missing values uses the SSRIH_grid mask between 1955 and 2015/2014. We compare the global land (except for Antarctica) or regional annual anomaly variation of SSR_20CR or SSR_CMIP6. The results show that SSR_20CR is significantly more consistent with the validation set than SSR_CMIP6.

Figure 6 Reconstruction capabilities of the AI model.

Figure 6a shows that the RMSE and CC of the SSR_20CR (0.247 W m${}^{-\mathrm{2}}/$0.970 W m⁻²) are smaller or larger than those of the SSR_CMIP6 (0.518 W m${}^{-\mathrm{2}}/$0.937 W m⁻²) with the original 20CR and CMIP6 dataset. The 20CR-AI model has a better reconstruction ability for SSR on the global land (except for Antarctica) scale. The RMSEs of the SSR_20CR (SSR_CMIP6) are 1.460 (2.413) W m⁻², 1.109 (1.829) W m⁻², 2.219 (2.596) W m⁻² and 1.286 (2.235) W m⁻² in North America, Europe, Asia and the Northern Hemisphere, whereas these values are 1.116 (1.766) W m⁻², 0.622 (1.602) W m⁻², 1.877 (1.839) W m⁻² and 0.772 (1.679) W m⁻² in South America, Africa, Australia and the Southern Hemisphere, respectively, concerning the original 20CR and CMIP6 dataset. In other words, the RMSEs of the SSR_20CR are smaller than those of the SSR_CMIP6 for the original 20CR and CMIP6 dataset except for Australia. In addition, the CCs of the SSR_20CR (SSR_CMIP6) are 0.958 (0.830) W m⁻², 0.958 (0.987) W m⁻², 0.886 (0.669) W m⁻², 0.930 (0.965) W m⁻², 0.938 (0.930) W m⁻², 0.943 (0.916) W m⁻², 0.936 (0.875) W m⁻² and 0.903 (0.822) W m⁻² in North America, Europe, Asia, the Northern Hemisphere, South America, Africa, Australia and the Southern Hemisphere, respectively, with respect to the original 20CR and CMIP6 dataset. That is, the CCs of the SSR_20CR are larger than those of the SSR_CMIP6 with the original 20CR and CMIP6 dataset except for Europe.

Based on the above comparison, the higher uncertainty for CMIP6 model output possibly biases the CMIP6-AI method. Thus, the accuracy of the SSR_20CR is higher than that of the SSR_CMIP6 at both global land (except for Antarctica) and regional scales. Therefore, we choose the reconstruction results of the 20CR-AI model as the final AI reconstruction dataset, and subsequent analysis in the following sections is only based on this dataset.

5.2 Comparison of the spatial and temporal variation characteristics

We investigate the long-term trends and spatial and temporal variation of the SSRIH_20CR, compare the differences between the SSRIH_20CR and SSRIH_grid, and suggest that the area and magnitude of the high and low centres of the SSRIH_20CR are the same as those of the SSRIH_grid. The results of the global land (except for Antarctica) reconstruction are consistent with dimming and brightening; the global dimming is primarily dominated by decreasing trends in Asia, Europe, Africa and North America, whereas Europe and North America are contributors to the increasing trends.

Figure 7 shows the spatial distribution of the SSRIH_grid and SSRIH_20CR for the 3 months (July 1960, July 1980 and July 2000). Figure S6 (see the Supplement) displays the spatial distribution of the annual SSRIH_grid and SSRIH_20CR from 1955 to 2018. Figure 7 also shows that the area and the magnitude of the high and low centres in the SSRIH_20CR are the same as in the SSRIH_grid. The SSRIH_20CR has mainly positive anomalies in Africa and the Eurasian continent in July 1960, especially in India and the Middle East. Afterwards, India showed a continuous and steady decline in SSR. This confirms the well-known phenomenon of global dimming over India (Wild et al., 2009; Soni et al., 2016, 2012; Padma Kumari et al., 2007; Kambezidis et al., 2012). In Australia, the SSRIH_20CR is dominated by negative anomalies in July 1980 and positive anomalies in July 1960 and July 2000. In Greenland, the SSRIH_20CR shows a large positive anomaly during 3 months. In northern Russia, there is a high value in July 2000. The reconstruction can better reflect the anomaly distribution of observation information, and the grid boxes with the missing values are infilled and reconstructed, which has high reliability.

Figure 7 Spatial distribution of the SSRIH_grid (a1–3) and SSRIH_20CR (b1–3) in typical months: 1–3 are July 1960, July 1980, and July 2000, respectively.

Figure 8 Global land (except for Antarctica) time series of the annual anomaly variations' SSR (relative to 1971–2000) before and after reconstruction.

Figure 8 illustrates the global land (except for Antarctica) annual anomaly variation and long-term trend of the SSRIH_20CR for the periods of 1955–2018, 1955–1991 and 1991–2018. Table S3 in the Supplement demonstrates the trends of global SSR change evaluation for various data sources on different scales. Also, we compare the differences between the SSRIH_20CR and SSRIH_grid. The minimum value of the SSRIH_20CR occurred in 1991 (−2.411 W m⁻²). The decreasing trend of the SSRIH_20CR from 1955 to 1991 (−1.276 ± 0.205 W m⁻² per decade) is slightly lower than that of the SSRIH_grid (−1.776 ± 0.230 W m⁻² per decade). After that, the SSRIH_20CR turns to an increasing trend of 0.697 ± 0.359 W m⁻² per decade from 1991 to 2018. This suggests that the difference between the SSRIH_20CR and SSRIH_grid may be caused by the results observed in limited data coverage (such as in Africa and North America) (Fig. 9). After homogenization and reconstruction, the trend (−1.276 W m⁻² per decade) from 1955 to 1991 corresponds to an overall reduction of −4.6 W m⁻² over the dimming period, while that (0.697 W m⁻² per decade) from 1991 to 2018 corresponds to an overall increase of 2 W m⁻² over the brightening period. This is in amazing agreement with the −4 W m⁻² for the dimming period and the 2 W m⁻² for the brightening period based on an overall surface energy budget assessment (Wild, 2012; see their Fig. 1). Also, similar conclusions (incomplete coverage of observational data leads to an underestimation of global warming trends) have been confirmed in global warming research (Gulev et al., 2021; Li et al., 2021).

Figure 9 Same as Fig. 8 but for regional annual anomaly variations.

Figure 9 demonstrates the long-term annual anomaly variations of the SSRIH_20CR in different regions and its results compared to the SSRIH_grid. Table S4 in the Supplement shows the evaluation in continental and hemispheric SSRIH_20CR change trends on different scales. The SSRIH_20CR shows a similar annual anomaly variation to the global land (except for Antarctica) average trend in North America and Asia, reaches a minimum in the late 1970s or early 1990s, and follows a moderate reversal. In Europe, the SSRIH_20CR shows a decrease (−2.180 ± 1.866 W m⁻² per decade) between 1963 and 1978 before turning to brightening (1.081 ± 0.312 W m⁻² per decade). In South America and Australia (Southern Hemisphere), the SSRIH_20CR shows no significant variation. In Africa, the SSRIH_20CR has a dimming trend (−1.506 ± 0.496 W m² per decade) from the 1950s to the 1990s, after which it remains levelled off (0.340 ± 0.998 W m⁻² per decade). The SSRIH_20CR shows a decreasing trend (−1.457 ± 0.246 W m⁻² per decade) until the 1990s in the Northern Hemisphere and a brightening (0.887 ± 0.415 W m⁻² per decade) afterwards. The annual average anomaly variations in regions and globally show that Asia, Africa, Europe and North America are the four contributors to the global dimming, while Europe and North America are two major contributors to the brightening. This is in general agreement with the results obtained by previous machine learning (Yuan et al., 2021). In addition, the discrepancy between the SSRIH_20CR and SSRIH_grid is more significant in low-coverage areas (right) than in high-coverage regions (left). It is particularly pronounced before 1980 and in South America. This suggests that the limited surface observations are not representative of the continental variation in SSR.

The sources of error in the observational dataset can be divided into three types. (1) Station errors are the uncertainties of individual station anomalies, including measurement errors (which are not the focus of the considerations in this paper) and errors due to homogenization. The errors due to homogenization adjustment are always approximately normally distributed (Jones et al., 2008; see their Fig. 5 and Fig. S9 in the Supplement) and therefore have limited impacts on the global average SSR change (Fig. S5a, b). (2) Sampling errors are the uncertainties in a grid box mean caused by estimating the mean from a small number of point values (Jones et al., 1997). (3) Bias error generally refers to systematic errors such as urbanization, which has not been discussed here. However, even the sum of the above errors is much smaller than the errors due to limited data coverage (Li et al., 2010; see their Fig. 5). So, the focus of this study is to eliminate this kind of error through the CNN reconstruction.

6 Data availability

Both the SSRIH_grid (the homogenized monthly gridded SSR data over 1923–2020) and the SSRIH_20CR (the monthly 20CR-AI model reconstructed SSR data for 1955–2018) are currently publicly available on the figshare website at https://doi.org/10.6084/m9.figshare.21625079.v1 (Jiao and Li, 2023). These datasets are also available at http://www.gwpu.net (last access: May 2023) for free.

7 Conclusions

In this study, we integrate global station observations based on the raw observational SSRs from GEBA and WRDC, combined with existing homogenized SSR datasets from other scholars. Also, we homogenize the globally distributed station data using the RHtestV4 software package. An improved CNN deep-learning algorithm is subsequently used to reconstruct the SSR anomalies. Thus, a reconstructed SSR anomaly dataset, SSRIH_20CR, is obtained based on training sets (20CRv3) for the years 1955–2018, with a resolution of 5^∘ × 2.5^∘. The main results are as follows.

1. The first integrated and homogenized global SSR monthly dataset is developed, which contains 944 stations in total and covers the longest periods (from the 1920s to recent years). A 5^∘ × 5^∘ grid box version of the monthly SSR anomaly dataset is derived.

2. This paper develops 5^∘ × 2.5^∘ full-coverage monthly land (except for Antarctica) SSR anomaly reconstructed datasets based on the above observations, using 20CRv3 to train the AI model. Comparative validations and evaluations show that the SSRIH_20CR provides a reliable benchmark for global SSR variations.

3. On average, the global annual SSR variations based on the SSRIH_grid are not significantly different, except that the increasing (brightening) trend after 1991 is a little smaller for the latter. The short-term brightening SSR in Europe from the 1970s to the 1980s disappears at the regional scale. At the same time, the brightening SSR after the 1990s in Asia slowed or was delayed.