Planetary albedo change exacerbates surface warming: a perspective from cloud-type changes

Persistent global warming is modulated by cloud changes, yet the specific contributions and associated controlling factors remain inadequately quantified. Using CERES radiation data with a surface energy-balance framework, we quantify the contribution of cloud radiative changes to decadal surface temperature trends over 2002–2023. Cloud changes exert a weak net effect on global mean warming due to near-cancellation between shortwave warming and longwave cooling, but strongly modulate the meridional pattern of zonal-mean warming. Specifically, clouds enhance warming in low- and mid-latitudes while mitigating warming at high latitudes. This pattern is associated with a systematic redistribution of cloud occurrence from low-/mid-level cloud types to high-level optically thin clouds, which reduce planetary albedo and weaken cloud longwave emission. These changes exhibit hemispheric difference. In 30–60° N, the region contributing most to global warming, the decline in the cloud-reflected solar radiation is mainly driven by decreased cloud fraction, linked to elevated sea surface temperatures, aerosol reductions, and mid-tropospheric drying. In 30–60° S, reduced cloud reflectivity resulting from decreased cloud optical thickness and increased liquid droplet radius dominates. However, at high latitudes in both hemispheres, increased mid-/high-clouds and enhanced cloud reflectivity, driven by enhanced moisture, upper-tropospheric static stability and increased cloud optical thickness, lead to greater reflected solar radiation and reduced downwelling longwave radiation, thereby attenuating local warming. Our results establish a direct observational link between cloud-type changes, planetary albedo decline, and contrasting warming across latitude zones, characterizing the cloud-radiative changes in recent climate change.

1 Introduction

Over the past two decades, the Earth is experiencing unprecedented warming, which exacerbates global heatwaves, wildfires, droughts, and glacier melt, and poses increasing threats to human health and life (Ripple et al., 2024; Gould et al., 2025). The 2023–2024 period experienced frequent record-breaking global mean surface temperatures (Li et al., 2024a; Cheng et al., 2025). While such record-breaking temperatures may be influenced by internal variability (Li et al., 2024a; Blanchard-Wrigglesworth et al., 2025; Goessling et al., 2025; Minobe et al., 2025), such as the El Niño-Southern Oscillation (ENSO), the magnitude of these extremes would be improbable without an underlying and sustained warming trend driven by earth's energy accumulation. This motivates a quantitative assessment of how observed radiative changes have accompanied and modulated the recent persistent global warming.

While the role of longwave radiation changes driven by factors such as greenhouse gases and aerosols in surface warming has been extensively studied, the contribution of shortwave radiation has often been overlooked. Planetary albedo (PA), which regulates the absorption and distribution of solar energy, plays a critical role in controlling surface temperature (T_s). A record-low PA has been identified as a key driver of both recent extreme temperature anomalies (Goessling et al., 2025) and marine heatwaves (Dong et al., 2025), while its accurate representation has proven crucial for improving the predictive skill of climate models (Blanchard-Wrigglesworth et al., 2025). Clouds, as the dominant regulator of PA (Stephens et al., 2015; Jian et al., 2018; Loeb et al., 2021; Li et al., 2024b, 2025b), influence the amount of solar radiation reaching the surface and thereby affect T_s. Liu et al. (2024) have highlighted that, after accounting for near-surface air temperature feedbacks, cloud radiative effects are the primary physical process driving monthly T_s variability over land. Under global warming, large-scale circulation adjustments and their associated dynamical and thermodynamic changes may trigger adjustments in cloud type, coverage, and albedo, thereby altering cloud radiative properties. Several model simulations have emphasized changes in cloud feedbacks under global warming. For example, Ceppi and Hartmann (2016) demonstrated that shortwave cloud radiative changes substantially intensify meridional temperature gradients and dominate the poleward expansion of large-scale circulation systems. Similarly, Zelinka et al. (2020) found that enhanced positive low-cloud feedbacks in CMIP6 models, driven by decreased extratropical low-cloud amount and optical depth, are the primary cause of their higher effective climate sensitivity. In addition, observations also reveal the extensive changes in cloud properties in recent years. Several studies have reported the reductions in low cloud cover at low-mid latitudes (Loeb et al., 2024; Goessling et al., 2025; Minobe et al., 2025). From a macroscopic perspective, Tselioudis et al. (2025) indicated that the contraction of storm zones drives the observed increase in global solar absorption. However, the contribution of these cloud changes to long-term T_s trends remains unquantified, and there is no observational evidence to support which specific cloud property changes have been most responsible.

To address these gaps, this study quantitatively assesses the radiative contribution of clouds to decadal T_s trends over 2002–2023 using satellite radiation data within a surface energy-balance framework. By isolating the cloud component of PA and decomposing it into cloud fraction and intrinsic cloud reflectivity for different cloud types, we separate their respective roles and further identify the underlying mechanisms. This observation-based assessment provides a benchmark for evaluating the representation of cloud radiative effects in climate models and offers new insights into the cloud type changes that influence recent warming trends.

2 Data and Methodology

2.1 CERES data

This study uses radiation and cloud data from July 2002 to February 2023 obtained from the Clouds and the Earth's Radiant Energy System (CERES) project. Two primary CERES products are employed: the Energy Balanced and Filled (EBAF) Edition 4.2.1 to quantify cloud radiative impacts on T_s trends and decompose PA, and the FluxbyCldTyp (FBCT) Edition 4.1 to analyze cloud-type changes.

We use monthly mean observed incoming solar radiation (ISR) and reflected shortwave radiation (RSR) at the top of atmosphere (TOA), alongside modeled surface upwelling and downwelling shortwave and longwave radiation fluxes at 1° resolution from EBAF (Wielicki et al., 1996; Kato et al., 2018; Loeb et al., 2018). This product provides diurnally complete averages by integrating geostationary satellite corrections (Loeb et al., 2018). Moreover, it ensures temporal continuity across satellite transitions and applies objective constraint algorithm to achieve energy balance closure based on ocean heat storage estimates (Loeb et al., 2009).

The FBCT Terra+Aqua Edition 4.1 product provides cloud-type partitioned radiative fluxes and cloud properties across pressure and optical depth bins (Sun et al., 2022), derived from CERES Single Scanner Footprint (SSF) Edition 4A and Moderate Resolution Imaging Spectroradiometer (MODIS) data. Both single-layer and multi-layer portions of upper-level clouds are classified as upper-level clouds. In this study, we employ 1° gridded monthly mean regional clear-sky shortwave fluxes, along with shortwave fluxes, cloud fraction (CF), cloud visible optical depth (CVOD), liquid cloud particle radius (r_liq) and ice cloud particle radius (r_ice) for different cloud types. Following the cloud type classification from International Satellite Cloud Climatology Project (ISCCP), we classify the clouds into nine types based on cloud top pressure and optical depth: cumulus (Cu), stratocumulus (Sc), stratus (St), altocumulus (Ac), altostratus (As), nimbostratus (Ns), cirrus (Ci), cirrostratus (Cs), cumulonimbus (Cb).

The shortwave cloud radiative forcing at the TOA (SWCRF_TOA) for the specific cloud type c is defined as:

$$\begin{array}{}\text{(1)}& {\mathrm{SWCRF}}_{\mathrm{TOA},c}={\mathrm{CF}}_{\mathrm{c}}\times \left({\mathrm{RSR}}_c-{\mathrm{RSR}}_{\mathrm{clr}}\right)\end{array}$$

where RSR_c is the upwelling overcast TOA flux for that cloud type, RSR_clr is the upwelling clear-sky TOA flux, and CF_c is the cloud fraction of that type. The term in parentheses represents the change in reflected shortwave radiation per unit cloud cover – i.e., the intrinsic shortwave cloud radiative effect of the cloud itself. We refer to this quantity as cloud reflectivity (CR):

$$\begin{array}{}\text{(2)}& {\mathrm{CR}}_c={\mathrm{RSR}}_c-{\mathrm{RSR}}_{\mathrm{clr}}\end{array}$$

CR is determined by the cloud's own microphysical properties, such as cloud thickness and the size of water droplets or ice crystals within the cloud. Multiplying CR by the corresponding cloud fraction yields the total shortwave cloud radiative forcing contributed by that cloud type.

The variable, SWCRF_TOA,c, is only used in decomposing the cloud component of RSR in Sect. 2.4. The classification thresholds and global mean climatology of CF and CR are provided in Fig. S1 in the Supplement.

2.2 Reanalysis data

In this study, monthly skin temperature data are obtained from the ECMWF Reanalysis v5 (ERA5). To further investigate the thermodynamic and dynamic drivers of cloud changes, several key monthly averaged meteorological variables from ERA5 are employed. Additionally, the aerosol optical depth (AOD) data from MERRA-2 are used to assess the potential influence of aerosols on cloud microphysics (Kaufman et al., 2005). The specific cloud controlling factors used are detailed in Sect. 2.5 and Table 1. All reanalysis datasets are regridded to a 1°×1° spatial resolution to match with the CERES data.

2.3 Radiative Contribution of Clouds to Surface Warming

Radiative forcing and climate feedbacks are often diagnosed using radiation fluxes at the TOA, which provide a useful framework for understanding global energy balance and temperature responses (Bony et al., 2006). However, regional T_s changes are strongly influenced by atmospheric and oceanic energy transports, so local temperature variations are not necessarily directly tied to TOA radiative anomalies (Liu et al., 2024). Therefore, examining the surface energy budget provides a more direct way to relate T_s changes to local radiative and turbulent energy fluxes (Andrews et al., 2009; Lu and Cai, 2009; Colman, 2015; Boeke and Taylor, 2018; Sejas et al., 2021), motivating the surface energy balance framework used in this study.

Following Lu and Cai (2009), we can quantify the contribution of clouds to T_s changes through both shortwave and longwave pathways. The monthly anomaly of T_s for each grid can be decomposed into contributions from different physical processes based on the surface energy balance framework (Lu and Cai, 2009):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Delta }{T}_{\mathrm{s}}& \approx \left\{-\left(\mathrm{\Delta }\mathit{\alpha }\right)\left(\stackrel{\mathrm{‾}}{{\mathrm{SW}}^{\downarrow }}+\mathrm{\Delta }{\mathrm{SW}}^{\downarrow }\right)+\mathrm{\Delta }{\mathrm{CRF}}_S+\left(\mathrm{1}-\stackrel{\mathrm{‾}}{\mathit{\alpha }}\right)\right.\\ & \left.\mathrm{\Delta }{\mathrm{SW}}_{\mathrm{clr}}^{\downarrow }+\mathrm{\Delta }{\mathrm{LW}}_{\mathrm{clr}}^{\downarrow }-\mathrm{\Delta }Q-\mathrm{\Delta }(\mathrm{SH}+\mathrm{LE})\right\}/\mathrm{4}\mathit{\sigma }{\stackrel{\mathrm{‾}}{T_s}}^{\mathrm{3}}\end{array}\end{array}$$

Each term on the right-hand side represents the contribution to T_s anomalies from distinct physical processes:

• (a)

  Surface albedo feedback (SAF): $-\left(\mathrm{\Delta }\mathit{\alpha }\right)\left(\stackrel{\mathrm{‾}}{{\mathrm{SW}}^{\downarrow }}+\mathrm{\Delta }{\mathrm{SW}}^{\downarrow }\right)$

• (b)

  Surface cloud radiative forcing:

  $$\begin{array}{}\text{(4)}& \mathrm{\Delta }{\mathrm{CRF}}_S=\left(\mathrm{1}-\stackrel{\mathrm{‾}}{\mathit{\alpha }}\right)\mathrm{\Delta }{\mathrm{SW}}_{\mathrm{cld}}^{\downarrow }+\mathrm{\Delta }{\mathrm{LW}}_{\mathrm{cld}}^{\downarrow }\end{array}$$
• (c)

  Non-SAF-induced clear-sky shortwave radiation: $\left(\mathrm{1}-\stackrel{\mathrm{‾}}{\mathit{\alpha }}\right)\mathrm{\Delta }{\mathrm{SW}}_{\mathrm{clr}}^{\downarrow }$

• (d)

  Clear-sky downward longwave radiation: $\mathrm{\Delta }{\mathrm{LW}}_{\mathrm{clr}}^{\downarrow }$

• (e)

  Surface heat storage (includes oceanic energy transport in oceans): −ΔQ

• (f)

  Combined sensible and latent heat fluxes: $-\mathrm{\Delta }(\mathrm{SH}+\mathrm{LE})$

Here, Δ denotes monthly anomalies relative to the climatological mean; overbars indicate monthly climatological means; α is surface albedo; SW (LW) denotes shortwave (longwave) radiation flux; downward arrows represent surface downward fluxes; SH (LE) represents sensible (latent) heat flux; and σ is the Stefan-Boltzmann constant. Subscripts clr and cld refer to clear-sky and cloudy conditions, respectively, with cloudy conditions calculated as the difference between all-sky and clear-sky conditions. Temperature units are in Kelvin (K).

Therefore, the monthly T_s anomalies caused by changes in surface cloud radiative forcing (ΔCRC, units: K) at each grid point is calculated as:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }\mathrm{CRC}=\left\{\left(\mathrm{1}-\stackrel{\mathrm{‾}}{\mathit{\alpha }}\right)\mathrm{\Delta }{\mathrm{SW}}_{\mathrm{cld}}^{\downarrow }+\mathrm{\Delta }{\mathrm{LW}}_{\mathrm{cld}}^{\downarrow }\right\}/\mathrm{4}\mathit{\sigma }{\stackrel{\mathrm{‾}}{T_s}}^{\mathrm{3}}\end{array}$$

The ΔCRC comprehensively captures both the shortwave radiative influences of clouds on T_s (ΔCRC_sw), which is not induced by surface albedo feedback, and the longwave radiative influences of clouds on T_s (ΔCRC_lw).

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathrm{\Delta }{\mathrm{CRC}}_{\mathrm{sw}}=\left(\mathrm{1}-\stackrel{\mathrm{‾}}{\mathit{\alpha }}\right)\mathrm{\Delta }{\mathrm{SW}}_{\mathrm{cld}}^{\downarrow }/\mathrm{4}\mathit{\sigma }{\stackrel{\mathrm{‾}}{T_s}}^{\mathrm{3}}\text{(7)}& {\displaystyle }\mathrm{\Delta }{\mathrm{CRC}}_{\mathrm{lw}}=\mathrm{\Delta }{\mathrm{LW}}_{\mathrm{cld}}^{\downarrow }/\mathrm{4}\mathit{\sigma }{\stackrel{\mathrm{‾}}{T_s}}^{\mathrm{3}}\end{array}$$

To analyze latitudinal patterns, we calculate area-weighted regional averages of ΔCRC within 30° latitude zones. The resulting time series for each latitudinal zone are then used to determine trends, representing the T_s trends caused by cloud radiative changes across different latitude zones (CRC_sw and CRC_lw).

The surface energy balance decomposition provides a useful diagnostic framework for quantifying the contributions of different energy flux components to T_s anomalies. However, its application over oceanic regions warrants careful interpretation. Over oceans, the large heat capacity of the ocean mixed layer and the influence of horizontal heat transport can result in a non-negligible local surface energy imbalance term (Donohoe and Battisti, 2013; Sejas et al., 2021). In addition, turbulent heat fluxes and ocean heat storage may covary with cloud radiative forcing. For instance, Liu et al. (2024) showed that over land, the cloud radiative effect is largely offset by its negative covariance with turbulent heat fluxes, as reduced cloud cover enhances surface heat loss. Therefore, the quantity ΔCRC derived from Eq. (5) should be interpreted as the direct radiative contribution associated with cloud changes under the linearized surface energy balance framework, rather than a complete attribution of T_s variability. Despite this limitation, the approach has been employed to diagnose the radiative effects of clouds on T_s at global scales, including over oceans (Izumi et al., 2015).

2.4 Decomposition of the cloud component of TOA RSR

In order to further clarify how clouds modulate solar radiation reaching the Earth's surface, we isolate the cloud component of TOA RSR, denoted as RSR_cloud, which represent the direct radiative effect of clouds on the TOA reflected shortwave radiation.

Following the simplified radiative model (Donohoe and Battisti, 2011; Stephens et al., 2015), the all-sky TOA RSR ($F_{\mathrm{TOA}}^{\uparrow }$) can be decomposed into the atmospheric and surface components (RSR_atm and RSR_surf).

$$\begin{array}{}\text{(8)}& F_{\mathrm{TOA}}^{\uparrow }={\mathrm{RSR}}_{\mathrm{surf}}+{\mathrm{RSR}}_{\mathrm{atm}}\end{array}$$

The simplified calculation equations are as follows (Li et al., 2025b):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{RSR}}_{\mathrm{atm}}={\displaystyle \frac{S\left(S{F}_{\mathrm{TOA}}^{\uparrow }-F_S^{\uparrow }{F}_S^{\downarrow }\right)}{S^{\mathrm{2}}-{F_S^{\uparrow }}^{\mathrm{2}}}}\text{(10)}& {\displaystyle }\begin{array}{rl}{\mathrm{RSR}}_{\mathrm{surf}}& =F_{\mathrm{TOA}}^{\uparrow }-{\mathrm{RSR}}_{\mathrm{atm}}=F_{\mathrm{TOA}}^{\uparrow }\\ & -{\displaystyle \frac{S\left(S{F}_{\mathrm{TOA}}^{\uparrow }-F_S^{\uparrow }{F}_S^{\downarrow }\right)}{S^{\mathrm{2}}-{F_S^{\uparrow }}^{\mathrm{2}}}}\end{array}\end{array}$$

where S is the incident solar radiation at TOA, $F_{\mathrm{TOA}}^{\uparrow }$ is the RSR at TOA, and $F_S^{\uparrow }$ and $F_S^{\downarrow }$ represent the upwelling and downwelling shortwave radiation fluxes at the surface, respectively. Similarly, the clear-sky RSR ($F_{\mathrm{TOA},\mathrm{clr}}^{\uparrow }$) can be decomposed into RSR_surfclr and RSR_atmclr.

$$\begin{array}{}\text{(11)}& F_{\mathrm{TOA},\mathrm{clr}}^{\uparrow }={\mathrm{RSR}}_{\mathrm{surfclr}}+{\mathrm{RSR}}_{\mathrm{atmclr}}\end{array}$$

The cloud component RSR_cloud is defined as the difference in the atmospheric component of RSR between all-sky and clear-sky conditions:

$$\begin{array}{}\text{(12)}& {\mathrm{RSR}}_{\mathrm{cloud}}={\mathrm{RSR}}_{\mathrm{atm}}-{\mathrm{RSR}}_{\mathrm{atmclr}}\end{array}$$

This quantity represents the solar radiation directly reflected by clouds themselves. Its magnitude and variability are primarily determined by the CF and CR. Notably, RSR_cloud is distinct from the traditional shortwave cloud radiative forcing at TOA (SWCRF_TOA), which also includes the effect of clouds masking the surface-reflected component. The SWCRF_TOA is defined as the difference in TOA RSR between all-sky and clear-sky conditions:

$$\begin{array}{}\text{(13)}& {\mathrm{SWCRF}}_{\mathrm{TOA}}=F_{\mathrm{TOA}}^{\uparrow }-F_{\mathrm{TOA},\mathrm{clr}}^{\uparrow }\end{array}$$

The relationship between the RSR_cloud and the SWCRF_TOA can be derived by substituting the decompositions:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathrm{SWCRF}}_{\mathrm{TOA}}& =\left({\mathrm{RSR}}_{\mathrm{surf}}+{\mathrm{RSR}}_{\mathrm{atm}}\right)-\left({\mathrm{RSR}}_{\mathrm{surfclr}}\right.\\ & \left.+{\mathrm{RSR}}_{\mathrm{atmclr}}\right)=\left({\mathrm{RSR}}_{\mathrm{atm}}-{\mathrm{RSR}}_{\mathrm{atmclr}}\right)\\ & +\left({\mathrm{RSR}}_{\mathrm{surf}}-{\mathrm{RSR}}_{\mathrm{surfclr}}\right)\\ & ={\mathrm{RSR}}_{\mathrm{cloud}}+F_{\mathrm{mask}}\end{array}\end{array}$$

where $F_{\mathrm{mask}}={\mathrm{RSR}}_{\mathrm{surf}}-{\mathrm{RSR}}_{\mathrm{surfclr}}$ represents the cloud masking effect on surface-reflected radiation – i.e., the reduction in TOA upwelling shortwave flux caused by clouds obscuring the surface.

Therefore, RSR_cloud can be equivalently expressed as the difference between the shortwave cloud radiative effect at TOA and this masking effect. Using the decomposition of SWCRF_TOA into contributions from nine cloud types (Eq. 1), we have:

$$\begin{array}{}\text{(15)}& {\mathrm{RSR}}_{\mathrm{cloud}}={\mathrm{SWCRF}}_{\mathrm{TOA}}-F_{\mathrm{mask}}=\sum _{c=\mathrm{1}}^{\mathrm{9}}\left({\mathrm{CF}}_c\times {\mathrm{CR}}_c\right)-F_{\mathrm{mask}}\end{array}$$

where CF_c and CR_c are the cloud fraction and intrinsic cloud radiative effect for cloud type c, respectively.

Thus, the monthly anomalies of grid RSR_cloud can be decomposed as follows:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}\mathrm{\Delta }{\mathrm{RSR}}_{\mathrm{cloud}}& =\sum _{c=\mathrm{1}}^{\mathrm{9}}\left(\mathrm{\Delta }{\mathrm{CF}}_c\times \stackrel{\mathrm{‾}}{{\mathrm{CR}}_c}+\mathrm{\Delta }{\mathrm{CR}}_c\times \stackrel{\mathrm{‾}}{{\mathrm{CF}}_c}\right.\\ & \left.+\mathrm{\Delta }{\mathrm{CF}}_c\times \mathrm{\Delta }{\mathrm{CR}}_c\right)-\mathrm{\Delta }{F}_{\mathrm{mask}}\end{array}\end{array}$$

The first two terms on the right-hand side represent the anomalies in RSR_cloud driven by changes in CF and CR, respectively, and the third term is their co-variation. The fourth term, −ΔF_mask, represents the anomalies in RSR_cloud arising from changes in the cloud masking effect (ΔF_mask). It is important to note that the cloud masking effect (ΔF_mask) is intrinsically tied to changes in CF and CR. Building on this decomposition, we can further investigate the sources of trends in RSR_cloud across different latitudinal zones. Note that in Fig. 4c, we plot the trends in −ΔF_mask term but not ΔF_mask to show its contribution to ΔRSR_cloud trends.

It is important to clarify that SWCRF_TOA used here differs from the surface shortwave cloud radiative effect (${\mathrm{SW}}_{\mathrm{cld}}^{\downarrow }$) employed in Sect. 2.3 to compute the cloud contribution to T_s trends; the latter is a surface flux quantity. SWCRF_TOA is used here solely for decomposing RSR_cloud and is not employed in subsequent analyses.

2.5 Quantifying Meteorological Contributions to Cloud Fraction Anomalies

To investigate the drivers of observed CF trends, we perform a stepwise multiple linear regression linking monthly CF anomalies for each cloud type to a set of relevant cloud-controlling variables. The cloud-controlling factors for low-level clouds have been extensively discussed in previous studies (Qu et al., 2015; Klein et al., 2017; Scott et al., 2020; Andersen et al., 2022; Naud et al., 2023; Li et al., 2025c; Naud et al., 2025), and those for high clouds have also been examined (Wilson Kemsley et al., 2024). In contrast, studies focusing on the controlling factors for middle-level clouds remain limited. Therefore, we select a set of candidate factors for middle clouds from dynamical, thermodynamical, and microphysical perspectives. Based on established physical understanding of cloud formation and dissipation, we select a set of cloud-controlling factors for each cloud height category (low, middle, and high clouds) as shown in Table 1. The detailed calculation methods for several variable selections are presented in Appendix A. All meteorological variables used in the calculations are derived from monthly mean ERA5 pressure-level data, while AOD is obtained from MERRA-2 monthly mean products.

Importantly, the observed CF trends and associated cloud-radiative changes may reflect a combination of aerosol-driven cloud adjustments (Bellouin et al., 2020), rapid adjustments to greenhouse-gas forcing and temperature-mediated cloud feedbacks (Zelinka et al., 2013; Sherwood et al., 2015), and internal climate variability (Zhou et al., 2016). Because each predictor may respond to more than one of these processes, the regression quantifies statistical relationships and trend contributions but does not uniquely attribute the observed trends to individual forcing, adjustment, internal-variability, or feedback components.

Table 1 Cloud controlling factors for low, middle, and high clouds.

Note: Factors labeled “This study” were not adopted directly from an established cloud-controlling-factor framework but were selected based on their expected dynamical or thermodynamic relevance to the corresponding cloud category. The accompanying references support the physical basis of the factor, although they may not use the identical variable definition or pressure level adopted here.

For each grid cell and cloud type, we first remove the seasonal cycles from CF and the meteorological predictors. To prevent common long-term trends from biasing the regression coefficients, linear trends are then removed from both the deseasonalized CF and predictor anomalies. The detrended predictor anomalies are standardized, whereas the detrended CF anomalies remain unstandardized. A stepwise multiple linear regression is then fitted as:

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }{\mathrm{CF}}_c^{\mathrm{DT}}\left(x,t\right)={\mathit{\beta }}_{\mathrm{0}}\left(x\right)+\sum _{k=\mathrm{1}}^K{\mathit{\beta }}_k\left(x\right)\cdot \mathrm{\Delta }{V}_k^{\mathrm{DT}}\left(x,t\right)\end{array}$$

where $\mathrm{\Delta }{\mathrm{CF}}_c^{\mathrm{DT}}$ is the detrended CF anomaly for cloud type c at location x and time t, $\mathrm{\Delta }{V}_k^{\mathrm{DT}}$ is the detrended and standardized anomaly of the kth meteorological variable, β_k are the regression coefficients, and β₀ is the intercept. The stepwise procedure is implemented with an entry significance level of p<0.05 and a removal significance level of p<0.10, ensuring that only variables with statistically meaningful explanatory power are retained. Because only the predictors are standardized, the resulting coefficients β_k represent the change in CF per unit standard deviation (SD) change in each predictor.

To estimate each factor's contribution to the observed CF trend, we apply the regression coefficient to the deseasonalized predictor anomaly before detrending, thereby retaining its observed long-term trend:

$$\begin{array}{}\text{(18)}& {\mathrm{Contrib}}_k\left(x,t\right)={\mathit{\beta }}_k\left(x\right)\cdot \mathrm{\Delta }{V}_k^{\mathrm{OBS}}\left(x,t\right)\end{array}$$

where $\mathrm{\Delta }{V}_k^{\mathrm{OBS}}$ denotes the deseasonalized predictor anomaly retaining its linear trend and scaled using the same standard deviation as in Eq. (17). For each of the six latitude zones, cloud type, and meteorological variable, we compute the area-weighted zonal mean contribution at each month. For each resulting zonal-mean time series, we then estimate the linear trend (% per decade) using ordinary least squares, representing the contribution of each meteorological factor to the long-term CF trend in that latitude zone.

2.6 Quantification of TOA longwave cloud masking

The conventional TOA longwave cloud radiative effect (LWCRE) is defined as the difference between clear-sky and all-sky outgoing longwave radiation.

$$\begin{array}{}\text{(19)}& {\mathrm{LWCRE}}_{\mathrm{TOA}}={\mathrm{OLR}}_{\mathrm{clr}}-{\mathrm{OLR}}_{\mathrm{all}}\end{array}$$

It can also be partitioned among the nine cloud types. For cloud type c, its contribution is calculated as

$$\begin{array}{}\text{(20)}& {\mathrm{LWCRE}}_{\mathrm{TOA},c}={\mathrm{CF}}_c\times \left({\mathrm{OLR}}_{\mathrm{clr}}-{\mathrm{OLR}}_c\right)\end{array}$$

where OLR_c is the overcast outgoing longwave radiation for cloud type c, OLR_clr is the corresponding clear-sky outgoing longwave radiation, and CF_c is the fractional coverage of that cloud type. Summing LWCRE_TOA,c across the nine cloud types gives the total LWCRE_TOA.

Changes in LWCRE do not arise exclusively from changes in cloud properties. Clouds modify the TOA radiative response to changes in non-cloud variables. For example, the radiative effects of changes in surface and atmospheric temperature, water vapor, greenhouse-gas concentrations, and ozone generally differ between clear-sky and cloudy conditions because clouds partly obscure longwave emission from the surface and lower atmosphere. Consequently, changes in these non-cloud variables can produce a change in LWCRE even if the clouds themselves remain unchanged. This contribution is referred to as the cloud-masking effect (Soden et al., 2004, 2008).

An LWCRE anomaly can therefore be expressed as

$$\begin{array}{}\text{(21)}& \mathrm{\Delta }{\mathrm{LWCRE}}_{\mathrm{TOA}}=\mathrm{\Delta }{\mathrm{LWCRE}}_{\mathrm{cloud}}+\mathrm{\Delta }{\mathrm{LWCRE}}_{\mathrm{mask}}\end{array}$$

where the longwave cloud-masking component is

$$\begin{array}{}\text{(22)}& \mathrm{\Delta }{\mathrm{LWCRE}}_{\mathrm{mask}}=\sum _x\left(\mathrm{\Delta }{\mathrm{OLR}}_x^{\mathrm{clear}}-\mathrm{\Delta }{\mathrm{OLR}}_x^{\mathrm{all}}\right)\end{array}$$

and x represents surface temperature (T_s), atmospheric temperature (T_a), water vapor (H₂O), well-mixed greenhouse gases (GHGs, including CO₂, CH₄, N₂O, CFCs, and HFCs), and ozone (O₃). We quantified these components using the publicly available ERA5 partial radiation perturbation (ERA5–PRP) dataset (Raghuraman et al., 2023b), following the methodology described by Raghuraman et al. (2023a). The dataset was generated using a two-sided PRP method with ERA5 monthly atmospheric fields, observed greenhouse-gas concentrations, and the offline RRTMGP radiative-transfer model. For each non-cloud variable, its cloud-masking contribution was calculated as the difference between its clear-sky and all-sky TOA longwave radiative perturbations.

As the ERA5–PRP dataset is available only through December 2020, this analysis covers the common period from July 2002 to December 2020, rather than the full July 2002–February 2023 period used in the main analysis. We use this overlapping period to demonstrate the magnitude and importance of longwave cloud masking rather than to extend the masking estimate over the full study period. The cloud-type decomposition is therefore interpreted as the observed LWCRE associated with each cloud category, which includes both cloud-property changes and cloud-masking contributions.

Figure 1 (a) ERA5 T_s trends from July 2002 to February 2023 and the T_s trends caused by changes in surface shortwave and longwave cloud radiative forcing (CRC_sw and CRC_lw) (methodology detailed in Sect. 2.3). (b) Global and latitude-zone averaged trends in the all-sky surface downward shortwave radiation fluxes (${\mathrm{SW}}_{\mathrm{all}}^{\downarrow }$)and the cloud component of TOA RSR (RSR_cloud), (c) CF, and (d) CR. Note that the latitudinal-mean trends in (a)–(c) are weighted by their area to account for their contributions to the global mean; the sum of the zonal trends approximately equals the global trend. Only trends significant at the 95 % confidence level are shown.

3 Results and Discussion

3.1 Surface Warming Trends and the Radiative Contribution from Clouds

Figure 1a shows the global and latitude-zone averaged T_s warming trend along with the T_s trends caused by changes in the surface shortwave and longwave cloud radiative forcing. Note that the latitudinal-mean values are area-weighted, and thus the sum of the zonal-mean trends constitutes the global trend. From 2002 to 2023, the global average T_s rose significantly by 0.22 K per decade. While the fastest warming occurred in the high-latitudes of the Northern Hemisphere (NH) (Fig. S2), the NH mid-latitudes contributes most to the global warming trend (Fig. 1a), primarily over Europe, the North Pacific, and the western Atlantic off the southeastern coast of North America (Fig. S2).

Previous research emphasized the role of longwave radiation in surface warming, particularly from greenhouse gases like water vapor. While clouds also contribute to the greenhouse effect by absorbing and re-emitting longwave radiation (Arking, 1991), cloud changes over the past two decades have induced a longwave cooling effect on T_s trends (Fig. 1a). This implies that the longwave effects alone are inadequate to accurately understand clouds' full radiative contribution to T_s trends. In fact, the radiative influence of clouds on the surface is twofold, involving longwave warming and shortwave cooling effects (Arking, 1991; Hartmann and Doelling, 1991; Wood, 2012). Changes in cloud properties also modify the amount of shortwave radiation reflected by clouds, thereby influencing the solar radiation reaching the surface. Although the surface longwave warming effect has weakened, the shortwave cooling effect has diminished even more markedly. Consequently, clouds have exerted a strong positive (warming) shortwave contribution that outweighs the negative longwave contribution, leading to a net enhancement of the global T_s trend, with notable latitudinal differences (Fig. 1a). Globally, the opposing shortwave and longwave cloud contributions largely cancel, resulting in a modest net cloud-induced warming of approximately 0.013 K per decade (derived from the sum of CRC_sw and CRC_lw in Fig. 1a). This suggests that non-cloud components – such as clear-sky greenhouse gas forcing – are the primary drivers of the global mean T_s trend.

However, the near-cancellation at the global scale belies the critical role of clouds in affecting the meridional variation of zonal-mean warming. The individual shortwave and longwave cloud effects represent large, opposing radiative perturbations, quantifying them separately is therefore essential for characterizing their contributions to the observed meridional warming pattern. At low and mid-latitudes, these effects compete: changes in CRC_sw enhance surface warming, while the CRC_lw partially offset it. Notably, in the Southern Hemisphere (SH) low latitudes, the CRC_sw trend even exceeds the total T_s warming trend, which is partially offset by a negative CRC_lwtrend. In the mid-latitudes of both hemispheres, the trends in the CRC_sw account for 57 % (SH) and 47 % (NH) of their T_s warming trend, underscoring its important role in the warming. In contrast, at high latitudes, the effects combine, as both cloud shortwave and longwave effects contribute to negative T_s trends, indicating that cloud changes partially mitigate polar warming (Kay and Gettelman, 2009; Alkama et al., 2020; Sledd and L'Ecuyer, 2021). Therefore, while the global net cloud impact on the T_s trend is small, clouds exert a meridionally contrasting influence that is fundamental to understanding the observed pattern of climate change.

To further investigate the drivers of cloud-induced T_s trends, Fig. 1b shows the global and latitude-zone averaged trends in surface downwelling shortwave radiation under all-sky condition (${\mathrm{SW}}_{\mathrm{all}}^{\downarrow }$). While the ${\mathrm{SW}}_{\mathrm{all}}^{\downarrow }$ directly affects surface energy balance, its changes are fundamentally governed by cloud-induced changes in RSR at the TOA. We therefore also examine the trends in cloud component of TOA RSR (RSR_cloud) (Fig. 1b) and its key drivers – CF (Fig. 1c) and CR (Fig. 1d). In 30–60° N and 0–30° S, the shortwave warming contribution from clouds (Fig.1a) is primarily due to significant increased solar radiation reaching the surface (Fig. 1b), resulting from reduced RSR_cloud (Fig. 1b), which is driven by both decreased CF (Fig. 1c) and reduced CR (Fig. 1d). This finding is consistent with Goessling et al. (2025), who highlighted that the recent decline in PA, particularly the record-low value in 2023, was predominantly driven by reduced low-level cloud cover over NH mid-latitudes and tropics, substantially amplifying global warming. In 0–30° N, although decreased CR leads to significantly less RSR_cloud, the solar radiation reaching the surface shows no significant increase (Fig. 1b), likely due to enhanced absorption by atmospheric water vapor (Loeb et al., 2021). Consequently, the shortwave cloud contribution to T_s trends is not significant here (Fig.  1a). Differently, in 30–60° S, despite no significant trend in ${\mathrm{SW}}_{\mathrm{all}}^{\downarrow }$ (Figs. 1b, S3a), clouds exhibit a substantially positive shortwave contribution to the T_s trend (Fig. 1a). This is because a decrease in RSR_cloud related to reduced CR contributes to an increase in solar radiation reaching the surface, which offsets the concurrent reduction in clear-sky surface downwelling solar radiation (Fig. S3b). This clear-sky reduction may reflect increased atmospheric attenuation associated with greater aerosol extinction and stronger water-vapor absorption (Haywood and Boucher, 2000; Wild, 2009; Trenberth et al., 2009; Li et al., 2024b). Therefore, clouds exert a positive shortwave contribution to T_s trends by counteracting the observed clear-sky reduction. Conversely, at high latitudes, increased CF and enhanced CR result in more solar radiation being reflected, reducing the solar radiation reaching the surface (Fig. 1b) and thereby partially mitigating surface warming.

Figure 2 Global and latitude-zone averaged trends in longwave radiation fluxes from July 2002 to February 2023. (a) All-sky surface downward longwave radiation (${\mathrm{LW}}_{\mathrm{surf},\mathrm{all}}^{\downarrow }$), (b) clear-sky surface downward longwave radiation (${\mathrm{LW}}_{\mathrm{surf},\mathrm{clr}}^{\downarrow }$), and (c) the difference between all-sky and clear-sky conditions (${\mathrm{LW}}_{\mathrm{surf},\mathrm{cld}}^{\downarrow }$ =${\mathrm{LW}}_{\mathrm{surf},\mathrm{all}}^{\downarrow }-{\mathrm{LW}}_{\mathrm{surf},\mathrm{clr}}^{\downarrow }$). (d) All-sky TOA outgoing longwave radiation (${\mathrm{LW}}_{\mathrm{TOA},\mathrm{all}}^{\uparrow }$), (e) clear-sky TOA outgoing longwave radiation (${\mathrm{LW}}_{\mathrm{TOA},\mathrm{clr}}^{\uparrow }$), (f) TOA longwave cloud radiative effect (TOA LWCRE = ${\mathrm{LW}}_{\mathrm{TOA},\mathrm{clr}}^{\uparrow }-{\mathrm{LW}}_{\mathrm{TOA},\mathrm{all}}^{\uparrow }$). Note that latitudinal-mean trends are area-weighted to account for their contributions to the global mean; the sum of the zonal trends approximately equals the global mean. Only trends significant at the 95 % confidence level are shown.

Concurrently, changes in cloud properties can also affect trends in the longwave cloud contribution to T_s trends (CRC_lw). Fig. 2 shows the trends in surface and TOA longwave radiation fluxes. In the 30–60° N and 0–30° S regions, although clear-sky components (e.g., water vapor and greenhouse gases) lead to an increasing trend in surface downward longwave radiation (Fig. 2b), the significant decrease in CF (Fig. 1c) causes the cloud-induced change in surface downward longwave radiation (i.e., the difference between all-sky and clear-sky; Fig. 2c) to exhibit negative trends. This indicates a weakened cloud greenhouse effect, allowing more surface heat to escape to space and thereby partially offsetting surface warming (Fig. 1a). Similarly, Liu et al. (2025) identified a widespread negative longwave cloud feedback over low- and mid-latitude continents, where declining low-level cloud cover reduces downwelling longwave radiation and partially mitigates surface warming. However, in the NH high latitudes, despite a significant increase in total CF, the cloud-induced change in surface downward longwave radiation shows a decreasing trend (Fig. 2c). One possible explanation is the observed redistribution of cloud vertical structure. Because cloud-base temperature generally decreases with increasing cloud-base height, lower and warmer clouds typically emit more longwave radiation toward the surface than higher and colder clouds (Intrieri et al., 2002; Shupe and Intrieri, 2004; Viúdez-Mora et al., 2015). Consequently, the concurrent decrease in low-level clouds and increase in high-level clouds may contribute to the reduction in cloud-related surface downwelling longwave radiation. At the TOA, the outgoing longwave radiation (OLR) trends (Fig. 2d–f) further support these interpretations. Over 30–60° N and 0–30° S, the reduction in CF allows more longwave radiation to escape to space, resulting in an increase in outgoing longwave radiation, corresponding to a negative TOA LWCRE trend. In contrast, over the NH high latitudes, the increase in high clouds (with colder cloud tops) reduces the OLR emitted to space. Overall, the changes in cloud amount and reflectivity are both relevant to the contrasting shortwave and longwave cloud contributions. The cloud-type analysis below evaluates these contributions more directly.

Figure 3 Latitude-zone and global averaged trends (July 2002 to February 2023) for CF and CR over (a–b) lands and (c–d) oceans. The CF trends for each zone are area-weighted to represent global contributions. Only trends significant at the 95 % confidence level are shown.

3.2 The Key Role of Cloud-type changes and Its Potential Mechanisms

Given the distinct radiative properties of clouds at different altitudes and optical depths (Hartmann et al., 1992), systematic changes in their distribution and characteristics significantly influence the cloud radiative effects and their contribution to T_s trends. The signal of these changes, however, is often masked within analyses of total cloud properties. To clarify this issue, we first examine the observed trends in CF and CR across nine cloud types (Fig. S4). These trends vary substantially between land and ocean (Fig. 3). We then quantify how these CF and CR changes affect the reflected solar radiation by clouds (RSR_cloud). Figure 4 displays the global and latitude-zone averaged RSR_cloud trends driven by changes in CF and CR across nine cloud types.

Figure 4 The trends in domain-averaged RSR_cloud caused by the changes in (a) CF and (b) CR across nine cloud types, and (c) the cloud masking effect (−ΔF_mask), at global and latitudinal scales, sorted in descending order of magnitude. Note that the latitudinal-mean trends are weighted by their area to account for their contributions to the global mean; the sum of the zonal trends approximately equals the global trend. Only trends significant at the 95 % confidence level are shown.

Globally, the decline in RSR_cloud is predominantly driven by the reduced CR, with the CF playing a secondary role (Fig. 4a, b). Specifically, the global RSR_cloud reduction trend of −0.22 W m⁻² per decade (the sum of contributions from all cloud types, including non-significant ones) is attributed to the cumulative decrease in CR of most cloud types (Fig. 4b), whereas the reduction of −0.14 W m⁻² per decade arises from compensating CF changes among different cloud types (Fig. 4a). The CF contribution reflects concurrent decreases in low-level (St and Cu), mid-level, and Cb clouds, which are partly offset by increases in optically thin high clouds (Cs and Ci). This redistribution among distinct cloud types and altitude categories is associated with broader changes in atmospheric circulation and thermodynamic conditions. As a result, the negative RSR_cloud contributions from the displaced low-level, mid-level, and Cb clouds are partially offset by the positive contributions from the increasing high-level optically thin clouds. This change would also reduce the downwelling longwave emission from clouds to the surface, thereby confirming the observed longwave cooling contribution of clouds to the T_s trend (Fig. 1a). The potential role of cloud masking effects will be discussed later in this section.

Figure 5 Decadal trends of standardized cloud controlling factors over the period 2002–2023 (units: per decade). Black dots indicate trends that are statistically significant at the 95 % confidence level (p<0.05). The factor list includes all variables considered for low-, mid-, and high-level cloud analyses.

Figure 6 Contributions of meteorological factors to CF trends (% per decade) for nine cloud types across six latitude zones, sorted in descending order of magnitude. Black dots indicate the observed significant CF trend, while red dots denote the sum of factor contributions. Latitudinal means are area-weighted to account for their contributions to the global mean. Only trends significant at the 95 % confidence level are highlighted.

To understand the drivers of the observed CF trends, we perform a stepwise multiple linear regression linking monthly CF anomalies to a set of cloud-controlling factors (Table 1). The temporal correlation coefficient between observed and regressed CF anomalies exceeds 0.6 for Cu, Sc, As, Ci, Cs, and Cb over mid- and low-latitude oceans, and reaches above 0.9 in the tropics (Fig. S5). For other cloud types, the correlation remains above 0.5 in regions where they are climatologically prevalent, indicating reasonable explanatory power. The contributions of each factor to the CF trend are calculated at the grid-cell level (by multiplying the regression coefficient by the local factor trend) and then area-weighted to zonal means. For clarity, we have only presented the zonal-mean coefficients and factor trends as shown in Figs. S6 and 5, respectively. Figure 6 shows the area-weighted contributions of individual factors to the CF trend across nine cloud types and six latitude zones. The close agreement between the observed CF trend and the sum of all CCF contributions indicates that the model captures a reasonable portion of the trends. However, the underlying mechanisms are more complex, and the statistical relationships should not be overinterpreted as a complete physical explanation, particularly for cloud vertical structure and microphysical processes (Wang et al., 2022; Zhao et al., 2026). Nonetheless, the close agreement in Fig. 6 suggests that the selected factors represent the dominant large-scale drivers, especially for cloud types that dominate the RSR_cloud decline.

Figure 7 Latitude-zone averaged trends in cloud properties (July 2002–February 2023) for (a) logarithm of cloud visible optical depth (log CVOD), (b) liquid cloud particle radius (r_liq), and (c) ice cloud particle radius (r_ice). Only trends significant at the 95 % confidence level are shown.

In contrast, the trends in CR are not explicitly regressed against meteorological factors, as CR is more directly linked to cloud microphysical properties. Figure 7 presents the zonal-mean trends in CVOD and particle radii (r_liq, r_ice), which serve as the primary controls on CR and provide complementary insights into the mechanisms underlying the RSR_cloud changes.

The most pronounced reductions in CF contributions occur in the NH mid-latitudes and the SH low latitudes (Fig. 4a). In the NH mid-latitudes, these decreases are most pronounced over the oceans (Fig. 3), mainly attributed to St, Sc, Ns, and Cb. The CF decline of St and Sc, which contribute most to the reduction in RSR_cloud, is largely statistically explained by rising T_s and decreasing AOD (Figs. 5, 6). The warming-induced suppression of low clouds is consistent with the positive cloud–temperature feedback mechanisms (Klein et al., 2017; Myers et al., 2018). Meanwhile, the tight coupling between CF reduction and AOD decline explicitly highlights the impact of aerosols on cloud lifetime (Albrecht, 1989). The reduction in aerosol loading may decrease cloud condensation nuclei (CCN) and cloud droplet number concentrations (CDNC), favoring larger droplets and more efficient collision–coalescence and precipitation, which may in turn reduce cloud persistence and fractional coverage (Twomey, 1974, 1977; Albrecht, 1989; Rosenfeld et al., 2019; Bellouin et al., 2020; Szopa et al., 2023). For Ns, the decreasing trend is mainly attributed to reduced AOD and increased mid-tropospheric dryness, as indicated by the positive VPD₅₀₀ trend. The higher VPD reflects lower relative humidity and less saturated conditions, making it harder for cloud droplets to form and persist, thereby inhibiting cloud growth. The decline in Cb is jointly driven by decreasing UTRH, and reduced AOD, while the increase in tropospheric thickness partially offsets this decline. Although the zonal-mean attribution suggests a substantial contribution from T_s (Fig. 6), this relationship is primarily evident over land (Fig. S7b), whereas the dominant decrease in Cb occurs over oceanic regions (Figs. 3, S7a). Therefore, increasing T_s does not explain the principal oceanic Cb reduction. Furthermore, the decreases in Ns and Cb are statistically consistent with thermodynamic and dynamical changes accompanying large-scale circulation adjustments under global warming. Model projections have suggested that in a warmer climate, baroclinic activity in the NH mid-latitudes may weaken due to Arctic amplification (Hadas et al., 2023). As a result, the frequency and intensity of mid-latitude cyclones are expected to diminish, thereby weakening large-scale ascent and directly reducing storm-associated cloud occurrence. Moreover, Tselioudis et al. (2025) also provide observational evidence of a continuous contraction of mid-latitude storm cloud regions over the past 24 years.

In the SH low latitudes, the largest contribution to the RSR_cloud decrease comes from the reduction in Cb, followed by Ns, As, and St. The reduction in Cb is primarily controlled by decreasing trends in upper-level ascent (as indicated by ω₃₀₀), reduced UTRH, and increasing trends in tropopause temperature. Similar to the NH mid-latitudes, the T_s signal in the zonal-mean attribution is dominated by land, while the observed Cb reduction is primarily oceanic. Moreover, although standardized ω₃₀₀ exhibits a negative zonal-mean trend over 0–30° S (Fig. 5), this zonal mean does not directly determine its contribution to the Cb trend because the contributions are calculated at each grid point before spatial averaging. The spatial distributions in Fig. S8 show that the negative ω₃₀₀–Cb regression coefficients over the western equatorial Pacific coincide with significant positive local trends in ω₃₀₀, indicating weakened ascent. This spatial correspondence produces a negative contribution of ω₃₀₀ to the Cb trend in this region (Fig. S8c), which dominates the negative zonal-mean contribution shown in Fig. 6. The local weakening of ascent over the western equatorial Pacific would tend to suppress convective development and vertical moisture transport, thereby limiting Cb formation. While recent studies have emphasized a strengthening trend in the Walker circulation over recent decades (Wu et al., 2021; Wills et al., 2022), our result should not be interpreted as a weakening of the Pacific Walker circulation as a whole. Instead, the local signal may reflect regional variability or a zonal reorganization within the broader Walker circulation. Taken together, the decline in Cb over 0–30° S is accompanied by reduced ascent and moisture supply and by upper-tropospheric warming and drying (Figs. 5 and 6), which may reflect a broader reorganization of tropical circulation. Their overall sign is consistent with a northward shift of the zonal-mean Intertropical Convergence Zone (ITCZ) and a corresponding adjustment of the ascending branch of the Hadley circulation (Schneider et al., 2014; Loeb et al., 2025; Guo et al., 2026; Shrestha et al., 2026), which would reduce deep convection and Cb occurrence south of the equator.

Notably, the net contribution of total CF changes to the RSR_cloud trend differs in sign between the two mid-latitude zones (Fig. 4). In 30–60° N, CF changes contribute approximately −0.12 W m⁻² per decade, reinforcing the decrease in RSR_cloud; this negative contribution arises mainly from decreases in St, Sc, Ns, and Cb. By contrast, in 30–60° S, the net CF contribution is positive (approximately 0.03 W m⁻² per decade) and therefore partly offsets the CR-driven decrease in RSR_cloud. This positive contribution can be attributed to an increase in more reflective Sc, which compensates for the RSR_cloud decrease caused by reduced Cu and mid-level clouds. As shown in Fig. 6, the increasing EIS trend and decreasing instability (M parameter) trend suppress convective Cu while favoring the formation of the more stable Sc.

The decrease in RSR_cloud contributed by reduced CR is primarily concentrated in the mid- to low-latitudes, particularly across the SH, with the Sc contributing most significantly (Fig. 4b). In both the SH low latitudes and NH mid- to low-latitudes, the diminished reflectivity of Sc can be attributed to significant decrease in its CVOD (Fig. 7a). Furthermore, a pronounced increase in the r_liq over mid- and low-latitude regions reduces the total scattering surface area of cloud droplets, resulting in less RSR (Fig. 7b). In addition, the decline in CR in the NH mid-latitudes coincides with a decreasing AOD trend (Fig. 5), which is potentially linked to controls on industrial and shipping sulfur emissions (Bai et al., 2020; Cao et al., 2023; Hodnebrog et al., 2024; Yuan et al., 2024; von Salzen et al., 2025). The aerosol decline may lead to decreased CCN, decreased CDNC, and larger effective droplet radii (Fig. 7b), thereby lowering cloud albedo (Twomey, 1974, 1977; Li et al., 2018, 2025a; Cao et al., 2023; Hodnebrog et al., 2024; Yuan et al., 2024). Similarly, Bai et al. (2020) have reported a persistent decline in both CDNC and AOD along the east coast of the United States, the west coast of Europe, and the east coast of China from 2003 to 2017. Moreover, von Salzen et al. (2025) have observed that the decrease in marine CR over the North Atlantic and Northeast Pacific from 2003 to 2022 can be attributed to the reduction in sulfur dioxide and other aerosol precursors. Notably, across SH low and mid-latitudes, a reduction in CR across nearly all cloud types has significantly contributed to the negative trend in RSR_cloud, which is consistent with their significant reduction in CVOD (Fig. 7a). This widespread darkening of clouds occurs despite increasing AOD in the SH mid-latitudes (Fig. 5), suggesting that the rising aerosol burden has not been efficiently translated into CCN (Cao et al., 2023). This is likely because the increased AOD is dominated by natural aerosols (e.g., sea salt) and above-cloud smoke, which are less effective as CCN or reside above the cloud layer, limiting their influence on boundary-layer clouds (Cao et al., 2023).

In contrast to the low- to mid-latitudes, increased CF and enhanced CR have led to a rise in RSR_cloud at high latitudes in both hemispheres, thereby mitigating regional warming (Fig. 4a, b). The increase in CF is primarily attributed to optically thick mid- and high-level clouds. Although the model performance is relatively limited at high latitudes (Fig. S5), the regression still captures part of the statistical relationship between CF anomalies and the selected environmental factors. They should, however, be interpreted as statistical indications rather than as evidence that the underlying physical processes are fully represented by the model. In the NH high latitudes, the positive trend in high clouds (Cs) is mainly linked to enhanced SUT and increased tropospheric thickness. For mid-level clouds (Ns), the increasing trend is primarily associated with enhanced mid-tropospheric moisture supply (RH₇₀₀) and strengthened warm advection (T_adv,700), which promote large-scale ascent and stratiform cloud formation. In the SH high latitudes, the increase in Ns is primarily driven by positive trends in AOD, RH₇₀₀, and CAPE, indicating the combined effects of increased CCN, enhanced moisture supply, and strengthened convective instability. The increase in high clouds (Cs and Cb) is associated with increasing AOD and SUT, along with decreasing T_trop. Additionally, the increase in St is linked to rising AOD and enhanced near-surface wind speed, which promote boundary-layer cloud formation through enhanced turbulent mixing. These statistically derived relationships are consistent with large-scale environmental changes in the polar regions. Enhanced atmospheric moisture (Ding et al., 2022; Patel and Kuttippurath, 2023) stems from intensified evaporation due to a wider ocean surface resulting from sea ice retreat (He et al., 2019) and increased poleward moisture transport from lower latitudes (Woods and Caballero, 2016; Kim et al., 2017), creating conditions more favourable for cloud development. This process is further amplified by the poleward shift of cyclone and storm tracks, which brings extensive cloud systems to the polar region (Tselioudis et al., 2024). Simultaneously, the observed cloud type changes, characterized by a decrease in high-emissivity low-level clouds and an increase in low-emissivity high-level clouds, provides the key to understanding the longwave cloud contribution to polar T_s trends. This change explains why the longwave cloud radiative effect has weakened, contributing negatively to the surface warming trend in both polar regions (Fig. 1a), despite the increased total CF in the NH and the stable total CF in the SH (Fig. 1c).

Moreover, the increase in CR at high-latitudes is mainly driven by low- and mid-level clouds in the NH and mid-level clouds in the SH. This enhancement is largely due to a significant increase in CVOD (Fig. 7a), likely resulting from greater moisture availability and elevated in-cloud liquid water content. Satellite observations have revealed that, over the past two decades, Arctic clouds have shown a tendency to shift from the ice to the liquid phase, with increased liquid water content enhancing cloud albedo (Lelli et al., 2023). More recently, Wang et al. (2026) found that declining mineral dust may amplify the ice-to-liquid replacement of high clouds over the NH mid- and high latitudes, increasing cloud optical depth and partially buffering warming through enhanced shortwave reflection. Additionally, the expansion of open ocean area enhances emissions of sea salt and biogenic aerosols, increasing the availability of cloud condensation and ice nuclei (Fig. 7c) (Schmale et al., 2021; Twohy et al., 2021; Lapere et al., 2023). This promotes higher CDNC and reduces droplet size (Fig. 7b, c), contributing to brighter and longer-lasting clouds.

It is noteworthy that the contribution of CF and CR changes to the trend in RSR_cloud also encompasses the masking effect of clouds on surface-reflected radiation (ΔF_mask). The total contribution of this effect to RSR_cloud trends (−ΔF_mask), defined as the difference of the TOA RSR surface component between clear-sky and all-sky conditions, is shown in Fig. 4c (see Sect. 2.4 for detailed explanation). This term (typically positive) exhibits a significant negative trend, which likely arises from two interrelated surface-cloud processes: first, a reduction in surface albedo due to declining snow/ice cover (Wu et al., 2020) or vegetation greening (Chen et al., 2019) lowers the clear-sky surface reflection (Li et al., 2024b), thereby diminishing the potential masking effect of clouds; second, a decrease in CF or a reduction in CR directly weakens the masking intensity by allowing more surface-reflected radiation to escape under cloudy conditions.

We further examined the cloud-type-resolved trends in TOA LWCRE (Fig. S9). Across the low and middle latitudes, nearly all statistically significant cloud-type contributions show decreasing LWCRE trends, except for a small positive Cu contribution in the NH midlatitudes. For several cloud types, particularly the low- and mid-level clouds and Cb over the NH midlatitudes and SH tropics, the negative LWCRE trends coincide with decreases in CF (Fig.  S4a). However, the correspondence is not one-to-one. For example, Cs occurrence increases globally, and both Cs and Sc increase over 30–60° S, whereas their associated LWCRE trends remain negative. These contrasting signs demonstrate that changes in cloud amount alone cannot explain the LWCRE trends. The concurrent decreases in cloud reflectivity for many cloud types, especially over the SH midlatitudes (Fig. S4b), suggest accompanying changes in cloud optical properties, although shortwave reflectivity does not directly constrain longwave emissivity. At a given cloud amount, LWCRE can also vary with cloud-top altitude and temperature, as well as cloud optical depth and emissivity (Zelinka et al., 2012; Raghuraman et al., 2024). Because these cloud-property contributions are not explicitly separated here, they are considered possible explanations rather than quantitatively attributed drivers. Moreover, conventional LWCRE includes the masking of changes in non-cloud properties by clouds (Soden et al., 2004, 2008; Raghuraman et al., 2023a). Using the ERA5–PRP dataset, we quantified the masking effects of surface temperature, atmospheric temperature, water vapor, well-mixed greenhouse gases, and ozone over the common period from July 2002 to December 2020 (Fig. S10). The global longwave cloud-masking trend is negative, with major negative contributions from water vapor and well-mixed greenhouse gases partly offset by the surface-temperature contribution, and accounts for approximately 32 % of the magnitude of the total negative LWCRE trend. Although this assessment ends 26 months earlier than the main analysis period, it demonstrates that longwave cloud masking is substantial.

4 Conclusions

Motivated by the rapid recent warming, ongoing climate crisis, and the need to understand cloud changes, this study quantifies the cloud radiative contribution to surface warming over 2002–2023 using CERES data and a surface energy balance framework. Globally, cloud-radiative changes cause shortwave warming, partly offset by longwave cooling. These contrasting effects are associated with a broad redistribution of cloud occurrence across vertical levels, characterized by decreases in low-/mid-level and Cb clouds, together with increases in high-level optically thin clouds (Cs, Ci), which is associated with large-scale circulation adjustments (Tselioudis et al., 2025), thermodynamic feedbacks (Klein et al., 2017), and microphysical processes (Li et al., 2018; Gui et al., 2026). This redistribution is consistent with the reduction in the cloud component of PA and the decrease in longwave radiation emitted by clouds towards surface.

Regionally, both cloud changes and the relative CF and CR contributions exhibit hemispheric differences. In 30–60° N, the RSR_cloud decline stems mainly from reduced CF (especially St, Sc, Ns, and Cb), which is consistent with rising SSTs, reduced AOD, and mid-tropospheric drying. In contrast, in 30–60° S, it is mainly driven by reduced CR, partially mitigated by increased Sc and Cs resulting from enhanced atmospheric stability. Climatologically, greater cloud reflection in the SH mid-latitudes compensates for the higher clear-sky albedo of the NH, helping maintain hemispheric PA symmetry (Bender et al., 2017; Datseris and Stevens, 2021; Jönsson and Bender, 2023; Hadas et al., 2023). During the study period, however, the larger reduction in RSR_cloud over the NH mid-latitudes is not balanced by the SH mid-latitude cloud changes, suggesting that this compensation pathway did not fully operate on the observed short timescale (Diamond et al., 2024; Loeb et al., 2025). Moreover, in 0–30° S, decreased CF (high-/mid-level clouds) and reduced CR (low-level clouds) contribute equally, associated with weakened ascent velocity, upper-tropospheric drying, and decreased CVOD with increased r_liq, respectively. Conversely, across both polar regions, clouds mitigate warming through both shortwave and longwave radiation changes, linked to increased CF and CR, particularly for mid/high clouds, resulting from increased water vapor, AOD, SUT, and enhanced CVOD.

The observed reductions in CF and CR are consistent with the long-term decline in global PA reported over recent decades (Loeb et al., 2021, 2024; Goessling et al., 2025; Hodnebrog et al., 2024; Li et al., 2024b, 2025b). Looking forward, this trend is projected to persist, with continued declines in RSR expected across the 21st century under all emission scenarios, particularly under high-emission pathways where the rate of decline accelerates (Li et al., 2025b). These findings highlight that the observed cloud changes in fraction and reflectivity not only contribute to current surface warming but also imply a sustained reduction in PA that will amplify future climate change. Moreover, given that subtropical low clouds dominate unforced interannual-to-decadal fluctuations in global energy imbalance (Miyamoto et al., 2026), their projected long-term decline under warming may not only amplify the forced energy uptake but also modulate the amplitude of future internal variability.

It should be also noted, however, that our findings are derived from passive satellite measurements, which have inherent limitations in distinguishing clouds from aerosols and over highly reflective surfaces (Yamanouchi et al., 1987; King et al., 1992; Mahesh et al., 2004), and cannot provide vertical cloud profile information. Future research based on active satellite observations is necessary to advance the understanding of cloud type changes.

Appendix A: Cloud-Controlling Factors and Variable Definitions

This appendix provides the definitions and calculation methods for the cloud-controlling factors summarized in Table 1.

Estimating Inversion Strength (EIS) is calculated according to the formula:

$$\begin{array}{}\text{(A1)}& \mathrm{EIS}=\mathrm{LTS}-{\mathrm{\Gamma }}_{\mathrm{m}}^{\mathrm{850}}\left(Z_{\mathrm{700}}-Z_{\mathrm{LCL}}\right)\end{array}$$

where lower-tropospheric stability (LTS) is defined as the potential temperature difference between 700 hPa and the surface, ${\mathrm{\Gamma }}_{\mathrm{m}}^{\mathrm{850}}$ represents the moist adiabatic lapse rate at 850 hPa, and Z₇₀₀ and Z_LCL denote the geopotential height at 700 hPa and the lifting condensation level relative to the surface, respectively. Following Wood and Bretherton (2006), the surface relative humidity is assumed to be 80 % to simplify the estimation of surface dewpoint temperature, with Z_LCL computed following Georgakakos and Bras (1984).

Surface temperature advection (ST_adv) is computed as:

$$\begin{array}{}\text{(A2)}& {\mathrm{ST}}_{\mathrm{adv}}=-{\displaystyle \frac{u}{R_{\mathrm{E}}\cos \mathit{\phi }}}{\displaystyle \frac{\partial \mathrm{ST}}{\partial \mathit{\lambda }}}-{\displaystyle \frac{v}{R_{\mathrm{E}}}}{\displaystyle \frac{\partial ST}{\partial \mathit{\varphi }}}\end{array}$$

where u and v denote eastward and northward 10 m winds, φ and λ represent latitude and longitude, R_E is the Earth's radius, and ST denotes the 2 m air temperature. Positive/negative ST_adv indicates warm/cold advection.

Vapor pressure deficit at 500 hPa (VPD₅₀₀) is defined as the difference between the saturation vapor pressure (SVP) and the actual vapor pressure (AVP) at 500 hPa (Murray, 1966; Buck, 1981):

$$\begin{array}{}\text{(A3)}& {\displaystyle }\mathrm{VPD}=\mathrm{SVP}-\mathrm{AVP}=\mathrm{SVP}\times \left(\mathrm{1}-\mathrm{RH}\right)\text{(A4)}& {\displaystyle }\mathrm{SVP}=\mathrm{0.61078}\times e^{\frac{a\times T_{\mathrm{a}}}{T_{\mathrm{a}}+b}}\end{array}$$

here RH is the relative humidity, T_a is the air temperature at 500 hPa. The constants a and b depend on the phase of water: for T_a at or above 0 °C (liquid water), a=17.27 and b=237.3; for T_a below 0 °C (ice), a=21.875 and b=265.5.

Moisture flux convergence at 700 hPa (MFC₇₀₀) is calculated as the negative divergence of the horizontal moisture flux:

$$\begin{array}{}\text{(A5)}& {\mathrm{MFC}}_{\mathrm{700}}=-\mathrm{\nabla }\cdot \left(q\mathit{V}\right)\end{array}$$

where q is specific humidity (kg kg⁻¹) and $\mathit{V}=(u,v)$ is the horizontal wind vector at 700 hPa. Convergence (positive MFC) indicates accumulation of moisture, which favors cloud formation.

Temperature advection at 700 hPa (T_adv,700) is given by:

$$\begin{array}{}\text{(A6)}& T_{\mathrm{adv},\mathrm{700}}=-{\displaystyle \frac{u_{\mathrm{700}}}{R_{\mathrm{E}}\cos \mathit{\phi }}}{\displaystyle \frac{\partial T_{\mathrm{700}}}{\partial \mathit{\lambda }}}-{\displaystyle \frac{v_{\mathrm{700}}}{R_{\mathrm{E}}}}{\displaystyle \frac{\partial T_{\mathrm{700}}}{\partial \mathit{\varphi }}}\end{array}$$

where u₇₀₀ and v₇₀₀ denote eastward and northward wind components at 700 hPa, φ and λ represent latitude and longitude, R_E is the Earth's radius, and T₇₀₀ is the air temperature at 700 hPa.

Upper-tropospheric static stability (SUT) is computed following the method described by Wilson Kemsley et al. (2024). Temperature and pressure profiles are first interpolated to 100 vertical levels using cubic spline interpolation from standard pressure levels to resolve fine vertical gradients. At each interpolated pressure level p, the static stability S_p is calculated as

$$\begin{array}{}\text{(A7)}& S_{\mathrm{p}}={\displaystyle \frac{R_{\mathrm{C}}}{c_{\mathrm{p}}}}{\displaystyle \frac{T_p}{p}}-{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}p}}\end{array}$$

where T_p is the temperature at pressure p, R_C is the gas constant, and c_p is the specific heat capacity at constant pressure. Finally, SUT is defined as the vertical average of S_p over the interpolated levels spanning from the tropopause pressure plus 50 hPa to the tropopause pressure plus 200 hPa. The tropopause pressure is identified on a monthly-mean basis using the standard World Meteorological Organization (WMO) definition (Reichler et al., 2003).

300 hPa zonal wind shear (ΔU300) is defined as the vertical gradient of the zonal wind between 300 and 925 hPa (Wilson Kemsley et al., 2024):

$$\begin{array}{}\text{(A8)}& \mathrm{\Delta }{U}_{\mathrm{300}}={\displaystyle \frac{U_{\mathrm{300}}-U_{\mathrm{925}}}{z_{\mathrm{300}}-z_{\mathrm{925}}}}\end{array}$$

where U is the easterly wind speeds and z the geopotential height.