1 Introduction

With the improvement in living standards, there is a growing desire for a better atmospheric environment and more clear blue skies. Traditionally, a sunny day with few clouds means blue sky day. However, with increasing anthropogenic emissions, the sky might be less or not blue even when the weather condition is “sunny (defined by cloud cover)” (Fast et al. 2006; Yang et al. 2011; Lin et al. 2012; Li et al. 2016a). On the other hand, fewer air pollutions do not directly bring in more blue skies, considering days with higher cloud cover or precipitation. The inconsistency between dropping PM and unchanged low visibility events (Ding et al. 2016) can be seen even during the COVID-19 lockdown period (Huang et al. 2020b; Huang et al. 2020c; Liu et al. 2020). All mentioned above try to convey the non-synergy correlation between air pollution and blue days when considering the unfavorable meteorological conditions dominated by large-scale synoptic circulation (Bei et al. 2020a). Thus, there are reasons to believe that “sunny” or “clean” may not be sufficient for a blue sky.

A series of heavy air-pollution events in the early 2010s in China have raised people’s strong desire for better air quality as well as more blue skies (Zhang et al. 2014; Zhang et al. 2015). The State Council of China has issued a series of strict anti-pollution control, such as the Air Pollution Prevention and Control Action Plan (China State Council 2013) and Blue-Sky Protection Campaign (China State Council 2018). However, there may be a discrepancy about whether clean actions have brought in more blue days by studies focusing on air quality or cloud cover (Zhang et al. 2019c; Bei et al. 2020b; Dai et al. 2020). Therefore, lacking direct index, it is still a challenge to understand the change of blue sky days in China and its underlying mechanism.

The first challenge is how to define a blue sky day. The occurrence of the blue sky depends not only on meteorological conditions but also on environmental conditions. First, a blue sky day should be a sunny day, which can be identified by low cloud cover and precipitation records. Second, a blue sky day should have clean air. Previous studies show that good air quality can be represented by high atmospheric visibility in the long term (Wu 2004; Ma et al. 2014; Zhang et al. 2014; Liu et al. 2017). Thus, our previous study (Wang et al. 2019) defined a blue sky day via the following criteria: the days with no rain, low cloud cover ≤75th percentile,, and the dry visibility in 14:00 ≥15 km.

However, the above definition does not tell the degree of blue of a blue sky day. In the real world, the gradations of the blue sky days are different from each other via presenting different RGB values (Fig. S1), which is dominated by the discrepancy between meteorology and the environment. A continuous deep blue days during the APEC (Asia-Pacific Economic Cooperation) conference in Beijing during December 2014 (Liu et al. 2016; Sun et al. 2016; Zhang et al. 2016; Xu et al. 2018) attract people’ attention around the world. Thus, classifying the blue days into different degrees will help us understand the change of blue days better.

The second challenge is how to assess the relative contributions of meteorological conditions and environmental factors to the change in the blue sky days. Observed evidence (Panmao et al. 1999; Alexander et al. 2006; Zhou et al. 2009; Hu et al. 2013; Easterling et al. 2019) shows that East Asia has experienced prominent changes in rainfall, temperature, and atmospheric circulation during recent decades. Meantime, with the development of industry, the emissions of aerosol have also increased significantly (Kato and Akimoto 2007; Meinshausen et al. 2009; Zhang and Wei 2014; Millar et al. 2017). Many studies show that the changes in meteorological and environmental conditions are responsible for the increase of air pollution events in China during recent decades (Kaiser 2002; Li et al. 2014; Dufour et al. 2015; Dai et al. 2020). For example, the weakened monsoon may increase the aerosol concentrations over eastern China (Zhu et al. 2012; Li et al. 2016b) and amplify persistent haze in Beijing (Li et al. 2016b; Cai et al. 2017). Nevertheless, we still do not know the role of meteorological and environmental factors in blue sky day change in China.

This study aims to creatively classify blue days into 3 grades, analyze their spatial-temporal variations during 1980–2018 in China, and reveal the effect of air pollution prevention plans. The rest of the paper is organized as follows. Section 2 describes the data and methods sued. Section 3 shows evidence for deep blue sky days decrease in China especially in winter and discusses the possible reasons for the change. Section 4 includes a conclusion and some discussions.

The classified blue days could be a direct index to evaluate the efficiency of the pollution prevention and control measures quantitatively in recent years. As life satisfaction and happiness are closely related to deep blue sky (Welsch 2006; Luechinger 2010; Li et al. 2014), this work will help to the construction of people’s livelihood. Meanwhile, it could provide a reference to studies on the relationship between weather and disease such as lung cancer and abortion (Pope et al. 2002; Zhang et al. 2019a; Zhang et al. 2019b). Focusing on days suitable for outdoor activities, it also serves the tourist industry and selection of important events. Finally, studying the mechanism behind the deep blue days will be a referee to host mega-events like Olympic Winter Games and implement policies on energy conservation and emission reduction.

2 Data and methods

2.1 Data

The daily meteorological observation data, including wind speed at 10 m, relative humidity, low cloud cover, and visibility values at 14:00 BJT from 1980 to 2018, were derived from the National Meteorological Information Center of the China Meteorological Administration. The monthly European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis-Interim atmospheric reanalysis data from January 1980 to December 2018 with a horizontal resolution of 0.75°×0.75° (Dee et al. 2011) are used in the analysis of climatic conditions for blue sky days. The Community Emissions Data System (CEDS) for historical emissions from 1980 to 2014 ((Hoesly et al. 2018) is used as an auxiliary to study the anthropogenic contribution, which is available from https://github.com/JGCRI/CEDS/.

2.2 Data quality controls

To ensure the accuracy and continuity of data used, a series of strict quality controls have been implemented following earlier work (Wang et al. 2019). In particular, as manual visibility measurement is broadly replaced by automated detection techniques around 2013–2015, special treatment is required for visibility records to keep data consistency. To this end, we transform the visibility during the post-2013 period according to the formula below (Organization 2008; Pei et al. 2018; Zhang et al. 2020a) to maintain data consistency:

$$\frac{VIS\ auto}{VIS\ manual} = \frac{\left(1/k\right)\times \ln \left(1/0.05\right)\ }{\left(1/k\right)\times \ln \left(1/0.02\right)}\approx 0.766$$

where k represents the extinction coefficient, and 0.05 and 0.02 denote the thresholds taken by the automatic and human observer, respectively. The visibilitymanual and visibilityauto are the human and instrument measurements of atmospheric visibility, respectively. In the analyses, the maximum of human observing visibility has been set to 30 km to reduce the discrepancy between the two observed methods. We also use the Standard Normal Homogeneity Test (SNHT) (Alexandersson and Moberg 1997) to verify whether stations exhibit considerable discontinuity in annual average visibility, by which only station observations passing the 99% confidence level test have been retained. Earlier studies have compared manual and automatic observations (Fan et al. 2017; Chang et al. 2020; Zhang et al. 2020a), and their results reported that although the threshold has been modified during 2013–2015, the continuity of the long-term data is not affected (Yin et al. 2017; Pei et al. 2018; Wang et al. 2020). Recent work demonstrated that this data adjustment can improve the continuity and reliability of visibility data (Zhang et al. 2020a). After quality control, data at 378 monitoring stations remain for further investigation (Fig. 1).

Fig. 1
figure 1

Location and altitudes of 378 meteorological stations in China

3 Methods

3.1 Classification of blue days

The classification procedures of blue days are as follows. Following our earlier work that takes the low cloud cover of the individual station into consideration, two thresholds are selected to represent the blue days at different levels. One is the national mean visibility (21.52 km) in all blue days from 1980 to 2018 in China (T1), and the other is the multi-year average value of low cloud cover in blue days at each station (T2). With the defined thresholds, we build two corresponding criteria to check: (1) whether the visibility is larger than T1; (2) whether the low cloud cover is less than T2. Based upon the above two criteria, Chinese blue days (CBDs) can be divided into three categories.

  1. (1)

    Deep blue: If one day meets both two prerequisites, it is recognized as the deep blue day (CBD3), representing the most desirable rate of blue days for people.

  2. (2)

    Medium blue: If it fits only one of them, it is termed as medium blue days (CBD2).

  3. (3)

    Light blue (CBD1): If it satisfies neither of the two prerequisites, it is ranked as the light blue days (CBD1), suggesting that it merely reaches the lowest level: the weakest blue days.

Case studies have been conducted to test the accuracy of the classification. Almost all well-known blue events, such as the APEC blue and Parade blue fall into the CBD3 category (not shown), indicating that our classification is reasonable.

4 Methods

A linear-trend analysis is used to reveal the variations of the time series of variables. Pearson correlation, and the standardized multivariate linear regression model (SMLR), and the Student’s t-test are performed to verify the relationships between CBD and the influencing factors. The SMLR model has been used in many studies to link predictand and predictors to quantify the relative contributions of each variable (Tai et al. 2010; Zhang et al. 2014; Zhai et al. 2019; Pei et al. 2020). By using the SMLR, the standardized regression coefficients (Beta) can be compared among different variables, which represents the percentage of change in CBD with unique variations of each variable if all the other variables are unchanged after removing impacts of units and sample size (Wilks et al. 2006). In this study, the square of Beta is used as the explained variance; the closer this statistic is to 100%, the greater the contributions to the classified CBD.

The Earth System Modeling Framework is used for interpolating data into the 1°×1° grid to increase statistical robustness (Jones 1999). Besides, the K index is employed to represent the stratification instability of moist air in the lower and middle troposphere. The equation is as follows (Zhang et al. 2007):

$$K=\left({T}_{850}-{T}_{500}\right)+{T}_{d850}-{\left(T-{T}_d\right)}_{700}.$$

where T and Td represent air temperature and dewpoint, and the subscript numbers show the geopotential height levels.

5 Results

5.1 Climatology and long-term trends of classified Chinese blue days in China

Figure 2a–c illustrate the spatial distribution of the average classified CBD in China during 1980–2018. At nationwide scale, the annual average CBD1-3 in China are 17.83, 82.01, and 62.00 days(d)/year(y), respectively. Basically, CBD1 exhibits a roughly geographic homogeneous pattern throughout the entire country, and the relatively low centers of CBD1 (below 5 days a year) are located in Yunnan, east of Tibet, and Hainan provinces. In contrast, Qinghai and most parts of Tibet share the maximum CBD1 with a value beyond 65 d/y. Compared to CBD1, the CBD2 tends to increase gradually from about 50d/y in southeast China to more than 80d/y in northwest China. As expected, almost all the low CBD2 centers reside in Southern China due to relatively higher low-cloud cover. The minimum CBD2 center, below 14 d/y, is located in Hebei province. CBD3 increases from southeast China to northwest China with the maximum number in Xinjiang and eastern Inner Mongolia (over 70 d/y), and the minimum number in western Qinghai and southern Hebei (Fig. 2c, over 13.46 d/y, and 26.85 d/y).

Fig. 2
figure 2

Spatial distribution of annual mean CBD under different grades in China during 1980–2018 (unit days), (a) CBD1, (b) CBD2, (c) CBD3, and (d) time series of the anomalous CBD in China from 1980 to 2018 (unit days)

The temporal evolutions of different CBD anomalies during 1980–2018 are shown in Fig. 2d. Upward trends are identified for CBD, CBD1, and CBD2 at a rate of 0.25, 0.34, and 0.39 days/year (d/y), respectively. Nevertheless, CBD3 is dominated by a declining trend exceeding −0.48 d/y, indicating a downward trend since 1980. Obvious increases of CBD1-2 can be found from 2014 to 2018, followed by an apparent decreasing CBD3 contemporaneous. This reveals the instant effect of emission controls on increasing CBD1-2 but an insufficient effect on increasing CBD3.

Figure 3 reveals the spatial distribution of trend coefficients with a horizontal resolution of 1°×1° in annual mean classified CBD from 1980 to 2018. It is obvious that CBD1 reveals an upward trend in almost 91.28% of grids (Fig. 3a). The most striking upward trend appears in northwestern China, with a value beyond 1 d/y. As for CBD2 (Fig. 3b), there is a similar dominant increasing trend over China (over 74.0%). Nearly all grids in the south of 38° N are dominated by an ascending trend, with a conspicuous center exceeding 0.95 d/y in southeastern China (SWC). In contrast, decreasing trends mostly occur in the north of 30° N, such as western Inner Mongolia, Qinghai, and some stations in northeast China. In contrast to CBD1 and CBD2, apart from some regions in western China and northeast China, CBD3 displays a coherent significant decreasing trend (more than 82.7% of, Fig. 3c). A noteworthy downward trend can be found in eastern China, with the minimum value below −1.18 d/y.

Fig. 3
figure 3

Observed distributions of trend coefficients (d/y) between 1980 and 2018 in China for (a) CBD1, (b) CBD2, (c) CBD3, and (d) CBD3 in winter; the white dots indicate the grid passing the 95% confidence level. The boxes overlaid indicate the six representative regions, Northeast China (NEC,), Northwest China (NWC), Yangtze-Huai River Valley (YHRV), Southeastern China (SWC), Pearl River Delta (PRD), and Beijing-Tianjin-Hebei (BTH)

Among the four seasons, the decrease of CBD3 is most prominent in winter. Fig. 3d shows the trend of wintertime CBD3 in China during 1980–2018, which is similar to the trend of annual CBD3. The regional mean winter CBD3 is decreasing at a rate of −0.11 d/y in China, with minimum centers located in northwestern China, southwestern China, and eastern China. These results are consistent with the increase in wintertime air pollution and an obvious increase in winter haze in China (Feng et al. 2020; Lu et al. 2020; Zhang et al. 2020b). In the following section, we focus on analyzing the reasons for the long-term change of wintertime CBD3 in China.

6 The reasons for wintertime CBD3 change in China

Four variables considered in our definition, including low cloud cover (LCC), visibility, relative humidity (RH), and non-precipitation day (non-REP) are taken into consideration. The non-REP is calculated by counting the monthly number of non-rainy and non-snow days in winter. Table 1 shows the trend coefficients of selected meteorological variables and their explained variances (Beta square values). According to the CBD3 definition, there should be a positive relationship between visibility and CBD3, and a negative one between LCC and CBD3. The RH plays a more complicated role by affecting the optical characteristics of hydrophilic aerosols (Peng et al. 2009; Petters and Kreidenweis 2013).

Table 1 The trend coefficients of winter CBD3 and the associated meteorological parameters over China and some sub-regions, and the Beta square values are given in the right column (units: %)

In China, the growth of non-REP and visibility should contribute to the increase of CBD3. Nevertheless, reduced RH contributes to the decrease of the dry visibility (below −0.04 km/y) and finally reduces the CBD3 along with increasing LCC. The changes of RH, non-REP, LCC, and visibility account for 17.39%, 9.73%, 69.22%, and 1.04% of the area-mean CBD3 change in China, respectively.

Moreover, we divided China into different sub-regions, including northwest China (NWC), southwest China (SWC), the Yangtze-Huai River Valley (YHRV), Beijing-Tianjin-Heibei (BTH), Pearl River Delta (PRD), and northeast China (NEC). The locations of these sub-regions are shown in Fig. 2d. In NWC, with the most significant reduction below −0.45 d/y, the decrease of CBD3 is dominated by reducing visibility and increasing LCC concurrently, with explained variances over 13.10% and 46.24%, respectively. In BTH, under unfavorable RH, LCC, and visibility conditions, the growth of non-precipitation days leads to a decreasing trend below −0.02 d/y, with a dominant role of RH beyond 95.84%. In YHRV, with unfavorable situations for all variables, an obvious decreasing trend (over −0.42 d/y) is detected. The changes in RH and LCC play a more important role in this region, accounting for 32.06% and 69.89% of the variations of CBD3, respectively. In SWC, where there is a unique downward trend (below −0.03%) of LCC, the decrease of RH and visibility contribute to a reduced CBD3.

In other regions with a faint increase (over 0.02 d/y) like PRD and NEC, variations of RH, visibility, and non-REP contribute to deep blue days, with a leading role (more than 30.58%) of non-REP and LCC (over 30%) in NEC. It should be mentioned that considering the close relationships among variations of meteorological factors, the combined multiple explained variances may be greater than 100% due to multicollinearity. Therefore, these results only give reference information; the underlined mechanism should be analyzed by numerical model in the future.

To reveal the long-term variations of CBD3, we analyze meteorological and environmental factors for the change of visibility, RH, and LCC. Previous studies indicate that meteorological elements could have a crucial impact dynamically and thermodynamically (Zhang et al. 2014; Liu et al. 2017; Huang et al. 2020a; Ma et al. 2020). Therefore, air temperature at 2 m (t2m), RH, WS, and K index are selected as major local meteorological conditions. The results are shown in Fig. 4a–d. A nationally downward trend of the K index in Qinghai, northwestern Sichuan, and southeastern Tibet can be found in Fig. 4a. A more stable atmospheric stratification in the lower and middle troposphere favors formations of LCC in relatively moist near-surface conditions (Wu et al. 2013; Horton et al. 2014). Meanwhile, the decrease of the K index should be related to the increase of t2m (Fig. 4b) and the decrease of RH (Fig. 4c). And it can partly explain the widespread increasing trend of CBD3, especially in coastal regions, consistent with previous studies (Cai et al. 2017; Li et al. 2018; Pei and Yan 2018). The wind speed also shows a significant decreasing trend in most parts of China (Fig. 4d). Previous studies show that decreasing wind speed can hinder the horizontal transportation of air pollutants (Zhang et al. 2014), which could partly explain the decrease of CBD3.

Fig. 4
figure 4

The spatial distributions of trend coefficients of (a) K index, (b) t2m (°C/year), (c) relative humidity (%/year), (d) wind speed (m/s∙year), (e) SO2, (f) NOx, (g) NMVOC, and (h) BC (units Tg/y) from 1980 to 2018 in winter. The white dots indicate the grid passing the 90% confidence level

Earlier works find that through scattering and absorbing function, aerosol extinction can directly deteriorate visibility, and sulfate ammonium, organics, and nitrate ammonium are the largest contributors to the light extinction coefficient in China (Huang et al. 2012; Yang et al. 2012). Meanwhile, anthropogenic aerosols could affect clouds mainly as the most common condensation nuclei (CNN) in environments by altering the microphysical properties of liquid and ice clouds (Bierkens et al. 2001). Additionally, increasing CNN can inhibit low cloud precipitation by weakening the mixed-phase stratocumulus (Yun et al. 2009; Qian et al. 2010). As a result, there exists an intimate linkage between anthropogenic aerosols and CBD3. Then, the trend of main precursors of fine particles (including SO2, NOx, and NMVOC) and BC (black carbon), which could affect visibility, are calculated. Only data from 1980 to 2014 are used due to data restriction.

Significant upward trends of SO2, NOx, NMVOC, and BC are shown in Fig. 4e–h displays, with a rate of 0.21, 0.21, 0.12, and 0.01 Tg/y during 1980 and 2014, respectively. The major emission-increasing zones are accompanied by obvious decreasing trends of CBD3, with the correlation coefficients of −0.43, −0.46, −0.47, and −0.51, all passing the 99% confidence level. The result reveals a stronger connection between CBD3 and human activity than meteorology at a national scale, and the BC is the primary emission affecting CBD3 in China.

YHRV, which is the most populous region with the highest emission (at rates over 0.84, 0.84, 0.46, and 0.04 Tg/y for SO2, NOx, NMVOC, and BC) in China, is found to be the most susceptible region. CBD3 in this region is significantly and negatively related to emissions, with the correlation coefficients all below −0.82. In NWC, SWC, and BTH, similar situations can be found except for the major role of BC in NWC and SO2 for another two, with a correlation coefficient of −0.39, −0.25, and −0.16, respectively. The weak connections might indicate the complicated feedback between aerosols and meteorology, which needs further study. On the contrary, weak correlations in PRD and NEC might indicate a more important role of meteorological conditions. We note that the results are probably affected by the limit of the analysis period. Further analysis is needed when the data are extended.

7 Conclusions and discussion

Although blue sky is a common weather phenomena and close to the needs of human, few studies have investigated it specially. In this study, we classify the types of blue sky according to meteorology condition (cloud cover and precipitation) and the concentrations of atmospheric aerosols (atmospheric visibility), and compute their long-term changes in China. It is found that deep blue sky days decrease significantly in China during 1980–2018, especially after 2013. NWC and YHRV share the most obvious downward trend on an annual scale or in winter.

Human activities at least partly contribute to the decrease of winter CBD3 in China. The variations of winter CBD3 are closely related to anthropogenic emissions of SO2, NOx, NMVOC, and BC with the correlation coefficients of −0.43, −0.46, −0.47, and −0.51, respectively. Meantime, the change of meteorological conditions is also responsible for the decrease of winter CBD3 in this period. The decrease in surface wind speed hinders air pollutant cleaning, and the increase in surface air temperature and decrease in relative humidity favor low cloud genesis, which is adverse to the occurrence of blue sky. These findings suggest that the occurrence of deep blue sky is affected by both the change of meteorological conditions and anthropogenic emissions.

This study only focus on China, but the definition of classified blue sky can extend to the global. Investigating how the blue sky change in the global may help us monitor the change of atmospheric environment of the earth, which deserve further study in the future.