Forest and bioenergy expansion amplifies climate warming by accelerating regional cloud loss

Land use and cover change (LUCC) can exacerbate cloud-mediated climate warming, yet long-term cloud responses to contrasting LUCC pathways remain poorly quantified. Here, we perform transient simulations based on the fully coupled Community Earth System Model (CESM) under the SSP1-2.6 scenario to quantify how long-term cloud trends respond to large-scale idealized afforestation (AF50, conversion of grassland or shrubland to forest), idealized bioenergy expansion (BE50, sugarcane expansion), and realistic afforestation (REAL). AF50 drives a pronounced vertical restructuring of cloud trends, steepening mid-level cloud decline and accelerating low-level cloud loss from −0.0153 ± 0.0030 % yr⁻¹ in the control simulation (CTL) to −0.0175 ± 0.0036 % yr⁻¹, corresponding to 1.14 times the CTL rate globally and 1.52 times over land. Concurrently, it amplifies the global increase in high-level clouds by a factor of 1.5 and accelerates total cloud cover loss over land. Bioenergy expansion induces a similar but weaker vertical shift, whereas realistic afforestation slows the loss rate of low- and mid-level clouds and suppresses regional warming. Interpretable machine learning identifies changes in relative-humidity trends as the highest-ranked predictor of LUCC-induced changes in low-level cloud trends, with relative importance values of 20.4 %, 20.4 %, and 15.3 % under AF50, BE50, and REAL, respectively. The simulated low-level cloud response is consistent with surface darkening reducing albedo, increasing sensible heat flux, deepening the planetary boundary layer, enhancing dry-air entrainment, and raising the lifting condensation level relative to the boundary-layer top, thereby suppressing low-level cloud formation. The cloud response is also modulated by remote atmospheric adjustment. The accelerated depletion of low-level clouds strengthens positive shortwave cloud radiative forcing, while the concurrent increase in high-level clouds enhances longwave warming, together amplifying regional warming hotspots. These results demonstrate that larger-scale forest expansion does not necessarily provide greater climate cooling benefits; the sign and magnitude of the outcome depend not only on the latitude band of land conversion but, more critically, on the type of land conversion.

1 Introduction

Large-scale atmospheric circulation systems on Earth directly facilitate cloud formation through thermodynamic uplift and cooling mechanisms (Boucher et al., 2013; Klein et al., 2017; Nomokonova et al., 2020). Clouds play a crucial role in modulating the radiative energy balance at the top of the atmosphere (TOA) and the Earth's surface, with their impacts indirectly linked to temperature variations (Schneider, 1972; Karl et al., 1993; Bretherton and Hartmann, 2009; Boehm and Thompson, 2023). For instance, low-level clouds, characterized by high albedo and significant optical depth, exert a substantial net cooling effect by reducing surface solar irradiance (Chen et al., 2000). In contrast, optically thin high-altitude clouds, such as cirrus, generally exert a weaker shortwave cooling effect because of their relatively low reflectivity, while their low cloud-top temperatures strongly reduce outgoing longwave radiation, thereby producing a net warming radiative effect (Liou and Ou, 1983; L'Ecuyer et al., 2019). The modulation of radiative balance and temperature by clouds is profoundly contingent upon their physical attributes, including type, thickness, liquid water path, and the dimensions of constituent particles (Albrecht, 1989; Stephens, 2005). The characteristics of cloud cover make it possible for them to play an amplifying role in climate warming (Ceppi and Nowack, 2021; Liu et al., 2025).

Clouds have been demonstrated to exert positive cloud feedback in nearly all climate models (Kiehl, 2007; Zelinka et al., 2013, 2017, 2020). Recent studies integrating satellite observations substantiate this assertion, revealing that cloud belts in equatorial and mid-latitude regions are contracting at a rate of 1.5 %–3 % per decade (Tselioudis et al., 2025). This subtle disparity between incident and reflected radiation, arising from variations in cloud fraction, exacerbates Earth's planetary energy imbalance and amplifies climate warming (Zhou et al., 2016; Li et al., 2024; Luo et al., 2024a; Mauritsen et al., 2025; Shi et al., 2025). Further evidence indicates that the diminution of low-level cloud cover in certain regions during the initial two decades of the 21st century has directly intensified anomalous regional temperature elevations, thereby facilitating the emergence of new climate warming hotspots (Tselioudis et al., 2024, 2025). Contemporaneous investigations elucidate that the recent widening of Earth's energy imbalance gap and record-breaking temperature highs are attributable to a marked decline in low-level cloud cover (Loeb et al., 2024; Goessling et al., 2025), with the underlying driver of cloud cover decrease being the escalation of greenhouse gas emissions from anthropogenic activities (Luo et al., 2024a). Among these, land use changes since the Industrial Revolution have contributed approximately one-third of the greenhouse warming through the emission of greenhouse gases (Houghton, 1999; Goldewijk and Ramankutty, 2004; Foley et al., 2005; van der Werf et al., 2009; Hong et al., 2021). Nevertheless, it remains unknown whether shifts in surface land cover types will render clouds less climatically beneficial and contribute to the formation of new warming hotspots.

Large-scale land use and cover change (LUCC) is altering circulation patterns in local and remote regions (Portmann et al., 2022). Forests, as a common form of land cover, are widely recognized for absorbing atmospheric carbon dioxide (CO₂) and mitigating climate warming (Peng et al., 2014). Changes in forests have been observed to exhibit significant interactions in heat fluxes and moisture between the surface and atmosphere from annual to seasonal scales; alterations in turbulence and water cycles will inevitably disrupt the atmospheric background for cloud formation and may reshape climate warming hotspots (Lian et al., 2022; Cheng et al., 2024). For instance, darker-colored forests absorb more solar radiation, leading to increased local radiation, whereby the absorbed additional energy is partitioned into sensible heat flux to heat the surface and the atmospheric temperature at the bottom of the boundary layer; this absorbed heat promotes the evaporation of surface and atmospheric moisture, thereby facilitating the development of shallow cumulus clouds (Ek and Holtslag, 2004). Compared to other vegetation, forests possess larger leaf area index (LAI) (Alkama and Cescatti, 2016), therefore, large-scale forest expansion enhances evapotranspiration by replacing croplands or grasslands, elevating atmospheric water vapor content and optimizing condensation environments, thereby increasing cloud cover and mitigating regional warming effects. Prior studies in Europe also indicate that forests overhead possess greater potential for cloud formation (Teuling et al., 2017). Another study on afforestation also demonstrates that the presence of forests augments low-level cloud amounts, and it is projected to yield positive climate benefits through such enhanced low-level cloud fraction; however, in regions with snow cover, this may result in diminished cloud amounts and engender converse climate effects (Duveiller et al., 2018, 2021; Breil et al., 2024). Recent studies have increasingly emphasized the critical role of cloud feedback in modulating the net climate impacts of land cover changes. For instance, Luo et al. (2024b) demonstrated that decreased cloud cover following deforestation partially offsets the cooling effect caused by increased surface albedo. Furthermore, Dror and Feingold. (2026) quantified how Amazon forest loss alters cloud properties and cloud albedo, exerting all-sky biophysical cooling feedback. These findings highlight the strong coupling between forest cover dynamics and cloud albedo responses, reinforcing the need to rigorously evaluate cloud-mediated radiative effects in land-use change scenarios.

Another distinctly divergent land use changes, namely the expansion of bioenergy crops, is less commonly implemented than forests owing to its comparatively diminished CO₂ sequestration capacity, yet it can engender sufficiently substantial climate mitigation extents through bio-geochemical and bio-geophysical processes (Georgescu et al., 2011; Wang et al., 2021; Cheng et al., 2024). However, conversely, the cultivation of bioenergy crops may precipitate unintended ramifications, such as the depletion of deep soil moisture and the amplification of water resource deficits (Righelato and Spracklen, 2007; Cheng et al., 2022; Li et al., 2023), with these adverse effects typically manifesting more pronouncedly at regional scales than at global scales (Hallgren et al., 2013). In contradistinction to the elevated evapotranspiration rates of forests, low-transpiration bioenergy crops (such as switchgrass) commonly exhibit lower LAI, a physical attribute that may diminish the efficiency of moisture conveyance between the atmosphere and terrestrial surfaces, thereby impeding enhancements in atmospheric humidity and consequently precipitating reductions in local cloud cover (Bala et al., 2007). Both forest and bioenergy expansions exert varying degrees of influence on temperature, water cycles, and clouds through radiative and non-radiative processes; such influences may render the atmosphere more humid or arid across ranges extending hundreds of kilometers, thereby modulating cloud formation and subsequently reshaping the configuration of warming hotspots via cloud radiative effect (Cerasoli et al., 2021; Wang et al., 2021). Although the impacts of LUCC on atmospheric circulation play a pivotal role in projecting future climate changes and have been extensively investigated (Feddema et al., 2005; Swann et al., 2012; Li et al., 2018), comprehension of its effects on long-term trends in cloud cover remains largely unknown, particularly as the expansion of bioenergy has not received commensurate attention to that afforded to forest expansion. Therefore, it is reasonable to hypothesize positing that LUCC affects the generation and trends of local or remote cloud cover through complex or straightforward mechanisms (Hua et al., 2023).

Large-scale LUCC has the potential to reshape future cloud cover trajectories via intricate bio-geophysical feedback (Luo et al., 2024b). To quantify these effects, we employed the Community Earth System Model (CESM) to conduct transient, fully coupled simulations spanning 2015 to 2100 (Danabasoglu et al., 2020). Our results demonstrate that while realistic afforestation may offer regional climate benefits, large-scale idealized forest and bioenergy expansion could exacerbate the decline of mid-to-low-level clouds while accelerating the increase of high-level clouds. These structural shifts in cloud cover amplify positive shortwave cloud radiative forcing, thereby intensifying the formation of regional warming hotspots. These results highlight the divergent climate impacts of land-cover transitions. While poorly constrained land-management scenarios risk undermining the cooling capacity of clouds and amplifying future climate warming, realistic afforestation pathways can effectively buffer this loss, thereby preserving essential biophysical cooling.

2 Data and Methods

2.1 Data

The ERA5 reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF), specifically the 2 m air temperature and cloud cover, were obtained from the Copernicus Climate Data Store (Muñoz-Sabater et al., 2021) to validate the CESM model simulations of near-surface temperature and cloud fraction. In addition, near-surface temperature fields from the National Centers for Environmental Prediction (NCEP) reanalysis were used to further evaluate the Community Earth System Model (CESM) model's temperature simulation performance and are available from https://psl.noaa.gov/ (last access: 1 September 2026) (Kalnay et al., 1996). The Global Precipitation Climatology Project (GPCP) dataset, used to validate the CESM model's performance in simulating recent precipitation, was obtained from https://www.psl.noaa.gov/data/gridded/data.gpcp.html (last access: 1 September 2026) (Adler et al., 2018; Huffman et al., 2023). An additional precipitation dataset is the CPC Merged Analysis of Precipitation (CMAP), which can be downloaded from https://psl.noaa.gov/data/gridded/data.cmap.html (last access: 1 September 2026) (Xie and Arkin, 1997). Finally, the core input data in this study is the land use data from Land-Use Harmonization 2 (LUH2) (Hurtt et al., 2020), downloaded from https://svn-ccsm-inputdata.cgd.ucar.edu/trunk/inputdata/ (last access: 1 September 2026).

2.2 Earth system simulation

We conducted the online simulations with the Community Earth System Model (CESM, version 2.1.3, https://www.cesm.ucar.edu/models/cesm2, last access: 1 September 2026) (Danabasoglu et al., 2020). As a state-of-the-art tool in climate prediction, this model accurately reproduces twentieth-century climate dynamics and surface sea temperatures while demonstrating high skill in projections of essential climate variables and atmospheric circulations (Danabasoglu et al., 2020; Fasullo, 2020). Similar to earlier versions (Lin and Rood, 1997), the atmospheric and land-surface component of CESM employs the Lin-Rood Finite-Volume (FV) dynamical core, though notable improvements have been introduced in atmospheric physics (Lin and Rood, 1997). Compared to its predecessor, the new generation FV mainly adjusted its cloud parameters (Zhu et al., 2022). CESM incorporates the Community Atmosphere Model Version 6 (CAM6) as its atmospheric component, offering a high-top option (with 70 pressure levels reaching about 130 km) and a standard option (with 32 pressure levels reaching about 40 km) (Danabasoglu et al., 2020). In this study, we utilize the standard low-top setup, which effectively simulates tropospheric atmospheric processes. CAM6 describes the atmospheric motion by the dynamics and sub-grid precipitation, cloud, and turbulence processes by the physics (Worley et al., 2011). Of particular relevance to our study is enhanced cloud-microphysics parameterization in CAM6: the Morrison–Gettelman 2 (MG2) scheme, which overcomes the prior constraint by explicitly predicting raindrop and snow-particle sizes (Morrison and Gettelman, 2008; Gettelman and Morrison, 2015). Accompanying enhancements feature the consolidated turbulence parameterization, incorporating a multivariate binormal probability density function (PDF) to collectively depict sub-grid fluctuations in temperature, humidity, and vertical velocity (Larson, 2017). At the same time, CESM classifies clouds across altitudes according to their cloud-top pressure, delineating low-level cloud from 700 hPa to the surface, middle-level cloud from 700–400 hPa, and high-level cloud above 400 hPa (Danabasoglu et al., 2020). Because high-level cloud is defined solely by cloud-top pressure, it is an aggregate diagnostic that may include physically distinct high-cloud types, including cirrus and the upper portions or anvils of deep convective clouds (Rossow and Schiffer, 1999; Bony et al., 2016). Accordingly, high-level cloud alone cannot identify the specific cloud type responsible for a simulated change or distinguish between surface-forced influences on deep convection and indirect large-scale dynamical or thermodynamic adjustment (Bony et al., 2016).

Before launching the simulations, we first evaluated the performance of CESM to reproduce historical climate change. The historical run employed the BHIST compset and was carried out online at 1° × 1° resolution. Monthly output in the period of 1979–2014 were used for model validation; to match the model grid, all observational datasets were conservatively interpolated to the CESM resolution. Annual and seasonal evaluations of temperature and cloud fraction show that CESM satisfactorily captures present-day climate variability (Figs. S1–S6 in the Supplement). More importantly, CESM reproduces the spatial patterns of cloud cover trends seen in the reanalysis data (Fig. S3), this is precisely the focus of this study. CESM also reproduces the observed global precipitation reasonably well, with a small global mean bias of only 0.2 mm d⁻¹ (Fig. S6). Then, the specific land-use forcing simulation design is implemented in the Community Land Model Version 5 (CLM5) within CESM. Relative to its predecessors, CLM5 incorporates advances in soil hydrology, crop modeling, fire feedback, and biogeochemical cycling, thereby enhancing the representation of anthropogenic land-use change within the climate system (Swenson and Lawrence, 2012; Li et al., 2017; Lawrence et al., 2019). The idealized simulation experiments in this study were all conducted based on Land-Use Harmonization 2 (LUH2) (Hurtt et al., 2020). CLM5 represents land-surface heterogeneity through a nested sub-grid hierarchy, allowing multiple land-unit to coexist within each grid cell (Oleson et al., 2013). Although CLM5 is not configured with a dynamic vegetation model, the time-varying landuse.timeseries dataset provided by the model allows vegetation's physiological processes (such as photosynthesis, transpiration, etc.) to be explicitly simulated (Lawrence et al., 2019; Cheng et al., 2022).

According to the Food and Agriculture Organization of the United Nations (FAO) definition of afforestation, it refers to the establishment of new forests on land that was not previously classified as forest through planting or deliberate seeding (FAO, 2012). This implies a fundamental change in land use types, for instance, from cropland to forest or from grassland to forest (FAO, 2012). Therefore, following the FAO definition, we constructed idealized land-use datasets for forest and bioenergy expansion simulations. In CLM5, there are eight tree plant functional types (PFTs), indexed from 1 to 8, and 64 crop functional types (CFTs). In CLM5, PFTs and CFTs coexist on nested sub-grid tiles that jointly mediate energy and moisture exchanges between the land surface and the atmosphere (Lawrence et al., 2019; Lombardozzi et al., 2020). We therefore implement forest and bioenergy expansion by altering their fractional areas at the sub-grid scale, consistent with previous studies (Dickinson and Henderson-Sellers, 1988; Swann et al., 2018). Specifically, our simulation experiments comprise a control run (CTL) in which land-cover fractions are held constant at 2015 levels, an idealized 50 % linear afforestation scenario (AF50), an idealized 50 % bioenergy crop expansion scenario (BE50), and a realistic forest expansion scenario (REAL) derived from previous studies (Roe et al., 2021; King et al., 2024; Weber et al., 2024), this land-use pathway exhibits a maximum 26 % forest-expansion pattern by 2100, representing an increase of 7.5 million km² relative to 2015 levels (Fig. S9), primarily distributed across tropical South America, Central Africa, and central North America (Figs. S7 and S8). Regarding the land area changes in the idealized afforestation scenario, approximately 29.3 million km² of grasslands and shrublands are converted to forests by 2100 relative to 2015 levels (Fig. S9), with a widespread global distribution (Figs. S7 and S8). Previous studies suggest that the global extent of grasslands and shrublands is sufficient to support land conversion of this magnitude (Di Vittorio et al., 2018); furthermore, this approach aligns more closely with the strict definition of afforestation established by the FAO. However, the REAL scenario does not strictly qualify as afforestation; it is more consistent with the definitions of reforestation or forest enhancement. Nevertheless, for analytical convenience, it is still referred to as an afforestation simulation in this study. All future-period simulations were conducted with the scientifically validated BSSP126 compset at 1° × 1° (f09_g17_gl4) horizontal resolution, while the ocean component was run at 0.5° × 0.5° resolution, the restart files for model initialization are derived from the National Center for Atmospheric Research's 500-year equilibrated pre-industrial control simulation (Eyring et al., 2016; Danabasoglu et al., 2020). In the idealized simulations, forest expansion is restricted to shrublands and grasslands, with the newly forested area distributed evenly among the eight tree PFTs, this practice of fundamentally converting land-cover types is also analogous to previous deforestation studies (Swann et al., 2018), wherein forests are replaced by grasslands or shrublands; bare ground is excluded, although some recent studies (Bastin et al., 2019; Wang et al., 2024; Liang et al., 2025) advocate afforestation on bare soil. Bioenergy expansion is represented by rained and irrigated sugarcane and occupies only existing bioenergy land, thereby encroaching neither on forests nor on natural grasslands. Bioenergy expansion is primarily driven by rainfed sugarcane cultivation. Relative to 2015, the cultivated area is projected to increase by 6.4 million km² by the end of the century (Fig. S9). Consistent with the Land Use Model Intercomparison Project (LUMIP) protocol (Lawrence et al., 2016), all experimental setups differ solely in their land-use forcing. Because atmospheric CO₂ is prescribed rather than prognostic, the simulated climate response originates entirely from biophysical processes induced by land cover change (e.g., changes in albedo and evapotranspiration). This configuration effectively excludes the offsetting effects of bio-geochemical feedback (i.e., carbon sequestration) on short-term climate response. The climate response is quantified as the difference between each sensitivity experiment and the control run (CTL), with statistical significance evaluated via a two-tailed Student's t test.

2.3 Attribution of LUCC-induced atmospheric changes to cloud trends

In this study, to disentangle the driving factors of cloud cover trends induced by different LUCC, we employed the explainable Extreme Gradient Boosting (XGBoost) model developed by Chen and Guestrin. (2016), whose proven capacity to resolve complex nonlinear interactions and deliver high predictive skill has made it a workhorse in both atmospheric and ecological modelling (Ren et al., 2024; Flora et al., 2024). After model training, we used SHapley Additive exPlanations (SHAP), based on game theory (Lundberg and Lee, 2017; Lundberg et al., 2020), to quantify the relative importance of individual predictors in the trained model and to characterize their relationships with model-predicted cloud-cover trends (Berdugo et al., 2022). SHAP values describe the contribution of each predictor to the model prediction relative to the model baseline and therefore should not be interpreted as a causal or process-based decomposition of the simulated physical cloud-cover change. For comparisons among predictors, global feature importance was calculated from the mean absolute SHAP value of each predictor and normalized by the summed mean absolute SHAP values of all predictors. The resulting percentages therefore represent normalized SHAP feature importance within the machine-learning model rather than the percentage of the physical cloud-cover change caused by a given variable. Specifically, we constructed 12 independent models for the cloud and climate-change trends driven by LUCC and performed 100 times Bayesian hyperparameter searches (Snoek et al., 2012) and the optimized parameters name and the obtained best values are listed in Table S1 in the Supplement. These predictors were selected to comprehensively represent the key physical processes controlling cloud formation and long-term trends under land use and cover change (Seneviratne et al., 2010; Stull, 1988; Wood, 2012; Santanello et al., 2018). They can be grouped into several physically meaningful categories: (1) surface energy and moisture fluxes, which directly reflect changes in surface energy partitioning induced by alterations in albedo and vegetation; (2) boundary layer structure, which captures the vertical extent and mixing intensity of the planetary boundary layer; (3) near-surface thermodynamic conditions, which govern saturation and cloud formation; (4) large-scale moisture availability and precipitation; and (5) atmospheric instability. In particular, sensible heat flux and planetary boundary layer height were included to represent critical upstream biophysical processes.

2.4 Calculation of cloud radiative effect

Cloud radiative effect (CRE) is commonly used to evaluate the cooling or warming influence of clouds on the Earth-atmosphere system (Goldblatt et al., 2021; Luo et al., 2024a), measured at the top-of-the atmosphere (TOA), in the atmosphere, or at the surface (Zhang et al., 2022), which consists of shortwave CRE (swCRE), longwave CRE (lwCRE) and netCRE (the sum of swCRE and lwCRE). The lwCRE was calculated as the difference between downward longwave radiation at the surface under cloudy conditions and that under clear sky (Luo et al., 2024a). The swCRE is determined by the difference between downwelling solar flux at surface (FSDS, W m⁻²) and clear-sky downwelling solar flux at surface (FSDSC, W m⁻²), as well as upward shortwave radiation at the surface (FSUS, W m⁻²) and clear-sky (FSUSC, W m⁻²), as indicated in Eq. (1) (Luo et al., 2024a).

$$\begin{array}{}\text{(1)}& \text{swCRE}=\text{FSDS}-\text{FSDSC}-\text{FSUS}+\text{FSUSC}\end{array}$$

In Eq. (1), FSUS and FSUSC are indirectly calculated via FSDS and FSDSC multiply albedo, as shown in Eq. (2). FSNS (W m⁻²) denotes the net shortwave radiation at the surface. The clear-sky albedo is also calculated using a similar equation (Wang et al., 2024).

$$\begin{array}{}\text{(2)}& \text{Albedo}=\mathrm{1}-\left({\displaystyle \frac{\text{FSNS}}{\text{FSDS}}}\right)\end{array}$$

2.5 Quantifying the impacts of land use change

We quantify the climate response to land-use change as the difference between the sensitivity simulations and the control run (CTL), represented by the symbol Δ. Using total cloud cover (CLDTOT, %) as an example, this effect is formally defined in Eq. (3).

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }\text{CLDTOT}=\text{AF50}-\text{CTL}\end{array}$$

2.6 Grid-based linear trend analysis

In this study, we conducted ordinary least squares (OLS) linear regression at every grid point to investigate the temporal variability patterns of atmospheric variables, encompassing cloud cover. The statistical significance of linear trend for each grid cell was derived using a Student's t test, as expressed in Eq. (4).

$$\begin{array}{}\text{(4)}& S={\displaystyle \frac{{\sum }_{t=\mathrm{1}}^n\left(x_{t,i,j}-{\stackrel{\mathrm{‾}}{x}}_{i,j}\right)\left(y_{t,i,j}-{\stackrel{\mathrm{‾}}{y}}_{i,j}\right)}{{\sum }_{t=\mathrm{1}}^n{\left(x_{t,i,j}-{\stackrel{\mathrm{‾}}{x}}_{i,j}\right)}^{\mathrm{2}}}}\end{array}$$

2.7 Quantification of the inequality of cloud cover trends

Originally developed as an economic metric for wealth distribution (Atkinson, 1970), the Gini coefficient (Gini index) is adapted in this study to assess the spatial disparity of hemispheric cloud cover trends, similar to its recent applications in studies of global precipitation changes and ocean productivity variations (Zhao et al., 2025; Lesk and Mankin, 2026). Equation (5) presents the continuous form of this index.

$$\begin{array}{}\text{(5)}& \text{Gini\hspace{0.17em}\hspace{0.17em}index}=\mathrm{1}-\mathrm{2}{\int }_{\mathrm{0}}^{\mathrm{1}}L\left(x\right)\mathrm{d}x\end{array}$$

In this formulation, the Gini index is derived by integrating the area beneath the Lorenz curve L(x), which quantifies the cumulative distribution of the target variable relative to the cumulative spatial area.

2.8 Surface energy balance decomposition of temperature

To better reveal the impacts of LUCC at different scales on surface thermodynamic temperature, this study employs the surface energy balance equation to dissect the various physical components of surface temperature changes. Through this decomposition method, the surface temperature changes caused by afforestation and bioenergy expansion can be attributed to multiple direct (local) or indirect (non-local) biophysical processes. Previous studies have shown that the linear approximation method via first-order Taylor expansion is sufficiently accurate for assessing surface temperature anomalies induced by large-scale land cover changes (Lee et al., 2011; Devaraju et al., 2018; Chen et al., 2020; Sui et al., 2025). According to the principle of energy conservation, the radiation energy absorbed by the surface must equal the heat released, as described in Eq. (6).

$$\begin{array}{}\text{(6)}& (\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}}){\text{SW}}_{\downarrow }+{\text{LW}}_{\downarrow }-\mathit{ϵ}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}=H+\mathit{\text{LE}}+G\end{array}$$

In Eq. (6), α_s represents surface albedo (dimensionless), and SW_↓ denotes downward shortwave radiation (W m⁻²); thus, $(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}}){\text{SW}}_{\downarrow }$ represents net shortwave radiation at the surface (W m⁻²). LW_↓ denotes downward longwave radiation at the surface (W m⁻²). ϵ is surface emissivity (dimensionless), a parameter used to describe how closely the Earth approximates a perfect blackbody, with typical values less than 1 (0.7 in this study). σ is the Stefan–Boltzmann constant (5.671 × 10⁻⁸ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$). H represents sensible heat flux (W m⁻²), LE represents latent heat flux (W m⁻²), T_s is surface temperature (°C), and G is surface heat flux (W m⁻², this value typically varies very little and can therefore be neglected in calculations). Rearranging and transforming the equation yields Eq. (7), where ρ is air density (1.29 kg m⁻³), C_p is specific heat capacity of air at constant pressure (1030 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), r_a is aerodynamic resistance (s m⁻¹), T_a is air temperature at 2 m (°C), and β represents the Bowen ratio (dimensionless).

$$\begin{array}{}\text{(7)}& \mathit{ϵ}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}=(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}}){\text{SW}}_{\downarrow }+{\text{LW}}_{\downarrow }-{\displaystyle \frac{(\mathrm{1}+\mathit{\beta })}{\mathit{\beta }}}{\displaystyle \frac{\mathit{\rho }{C}_p}{r_{\mathrm{a}}}}(T_{\mathrm{s}}-T_{\mathrm{a}})\end{array}$$

Linearizing Eq. (7) yields the difference in surface temperature changes ΔT_s between each sensitivity experiment and the CTL experiment. After linear decomposition, the equation allows the surface temperature change ΔT_s to be explicitly expressed as a function of albedo (α_s), turbulent processes (β), surface roughness or aerodynamic resistance (r_a), cloud radiation effects (LW_↓ and SW_↓), and atmospheric feedback (T_a), as shown on the right-hand side (RHS) of Eq. (8). The first three physical terms on the RHS represent the local changes induced by LUCC, while the latter three terms indicate non-local responses (Devaraju et al., 2018).

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{\Delta }{T}_{\mathrm{s}}=& \left({\displaystyle \frac{-{\mathit{\lambda }}_{\mathrm{0}}{\text{SW}}_{\downarrow }}{\mathrm{1}+f}}\right)\mathrm{\Delta }\mathit{\alpha }+\left[{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{0}}H}{(\mathrm{1}+f){\mathit{\beta }}^{\mathrm{2}}}}\right]\mathrm{\Delta }\mathit{\beta }\\ & +\left[{\displaystyle \frac{fH}{(\mathrm{1}+f)\mathit{\rho }{C}_p}}\right]\mathrm{\Delta }{r}_{\mathrm{a}}+\left[{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{0}}(\mathrm{1}-\mathit{\alpha })}{\mathrm{1}+f}}\right]\mathrm{\Delta }{\text{SW}}_{\downarrow }\\ & +\left({\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{0}}}{\mathrm{1}+f}}\right)\mathrm{\Delta }{\text{LW}}_{\downarrow }+\left({\displaystyle \frac{f}{\mathrm{1}+f}}\right)\mathrm{\Delta }{T}_{\mathrm{a}}\end{array}\end{array}$$

In Eq. (8), λ₀ and f represent the temperature sensitivity coefficient and the energy redistribution coefficient, respectively, as defined by Eqs. (9) and (10). These two parameters, respectively describe the response of black-body radiation to temperature and how turbulent fluxes regulate surface temperature (Dickinson, 1983; Lee et al., 2011).

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\sigma }\mathit{ϵ}{T}_{\mathrm{s}}^{\mathrm{3}}}}\text{(10)}& {\displaystyle }& {\displaystyle }f={\displaystyle \frac{\mathit{\rho }{C}_p}{\mathrm{4}\mathit{\sigma }\mathit{ϵ}{T}_{\mathrm{s}}^{\mathrm{3}}{r}_{\mathrm{a}}}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}\right)\end{array}$$

Finaly, we calculate the annual trend of each term using Eq. (4). In this way, the long-term trends of each physical term induced by LUCC forcing can be investigated.

2.9 Attribution of non-local temperature responses via horizontal temperature advection

Temperature advection is a primary driver of non-local temperature responses to large-scale LUCC. Changes in surface energy fluxes trigger boundary-layer dynamics and mesoscale circulations that transport thermal local perturbations through wind fields to remote regions (Pielke et al., 1998; Fig. S10). Following previous studies (Kan et al., 2025), we evaluate horizontal temperature advection (T_adv) at 850 hPa, the critical interface between the boundary layer and the free troposphere. This is defined in Eq. (11), where u and v are the horizontal winds at this standard isobaric level, and $\frac{\mathit{\vartheta }T}{\mathit{\vartheta }x}$, $\frac{\mathit{\vartheta }T}{\mathit{\vartheta }y}$ represent the zonal and meridional temperature gradients, respectively. Therefore, temperature advection can be decomposed into zonal ($T_{\text{adv}}^u$) and meridional ($T_{\text{adv}}^v$) components. Finally, the difference of temperature advection (ΔT_adv) between the sensitivity and control experiments elucidates the atmospheric dynamic mechanisms underpinning LUCC-induced teleconnections.

$$\begin{array}{}\text{(11)}& T_{\text{adv}}=T_{\text{adv}}^u+T_{\text{adv}}^v=\left(-u{\displaystyle \frac{\mathit{\vartheta }T}{\mathit{\vartheta }x}}\right)+\left(-v{\displaystyle \frac{\mathit{\vartheta }T}{\mathit{\vartheta }y}}\right)\end{array}$$

2.10 Diagnosing planetary wave activity using Plumb fluxes

To capture the dominant propagation characteristics of Rossby waves in the middle-to-upper troposphere, we evaluated the wave activity at the 500 hPa isobaric level, consistent with previous studies (Harada and Nakagawa, 2017; Kan et al., 2025). This altitude serves as a canonical diagnostic layer for vigorous planetary wave activity, particularly for quasi-stationary and stationary Rossby waves. Located near the mid-to-upper tropospheric jet core, wave activity fluxes at 500 hPa clearly manifest the mid-to-upper atmospheric response to surface perturbations, while mitigating the interference of near-surface friction and turbulent dissipation. Consequently, the wave activity flux (WAF) constitutes an essential tool in atmospheric dynamics for diagnosing Rossby wave energy transport, effectively tracing the propagation pathways, directions, and intensities of planetary waves. In this study, we utilized the two-dimensional WAF formulation (Plumb, 1985). Because this flux is parallel to the Rossby wave group velocity, it intuitively delineates the two-dimensional propagation of planetary wave energy, as detailed in Eq. (12).

$$\begin{array}{}\text{(12)}& F_{\mathrm{s}}={\displaystyle \frac{p}{p_{\mathrm{0}}}}\cos (\varnothing )\times \left[\begin{array}{l}{u}^{\prime \mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}\mathrm{\Omega }a\sin \mathrm{2}\varnothing }\frac{\partial \left(v^{\prime }{\mathrm{\Phi }}^{\prime }\right)}{\partial \mathit{\lambda }}\\ -u^{\prime }{v}^{\prime }+\frac{\mathrm{1}}{\mathrm{2}\mathrm{\Omega }a\sin \mathrm{2}\varnothing }\frac{\partial \left(u^{\prime }{\mathrm{\Phi }}^{\prime }\right)}{\partial \mathit{\lambda }}\end{array}\right]\end{array}$$

Here, the prime denote the deviation from zonal mean; ∅, λ, Φ denotes latitude, longitude, and geopotential height, respectively; 2Ωsin (∅) represents the Coriolis parameter, where a is the Earth's radius and Ω is the Earth's rotation rate. p₀ denotes the near-surface reference pressure (1000 hPa). The u and v denote the zonal and meridional winds, respectively. This set of notations and parameter configurations is consistent with the two-dimensional wave activity flux formulation has been validated in numerous climate modeling studies (Huang and Nakamura, 2016; Kan et al., 2025).

3 Results

3.1 Long-term response of cloud cover to LUCC forcings

Our simulations demonstrate that distinct land-use forcing significantly altered the global long-term trends in total cloud cover (CLDTOT) and low-level clouds (CLDLOW) (Fig. 1). While CLDTOT exhibited a moderate negative trend in both CTL (−0.0083 ± 0.0020 % yr⁻¹, R = −0.6547) and the bioenergy expansion experiment (BE50, −0.0083 ± 0.0020 % yr⁻¹, R = −0.6363), the idealized (AF50) and realistic afforestation (REAL) scenarios modestly attenuated this decline to −0.0072 ± 0.0025 and −0.0069 ± 0.0021 % yr⁻¹, respectively. This suggests that the enhanced evapotranspiration and increased surface roughness associated with afforestation, coupled with a compensatory increase in cloud cover at other cloud-top heights, may provide a buffering effect against the loss of total cloud cover (Duveiller et al., 2021; Xu et al., 2022). Secondly, this global attenuation is largely driven by a weakened cloud loss over the oceans, where the CLDTOT decline rate of −0.0073 ± 0.0019 % yr⁻¹ in the CTL simulation was mitigated to −0.0054 ± 0.0025 % yr⁻¹ in AF50 (Fig. S11a and b). Notably, however, idealized afforestation exacerbated the loss of CLDTOT over land, accelerating the decline rate by a factor of 1.05 relative to the CTL run. CLDLOW exhibited the most pronounced and contrasting response. While idealized afforestation accelerated the decline of CLDLOW by a factor of 1.14 relative to the CTL simulation, this effect was notably amplified over land, where the loss rate increased by a factor of 1.52 (Fig. S11e and f). This suggests that the extensive conversion of grasslands or shrublands into forests may directly reduce low-level clouds by altering surface evapotranspiration and promoting deeper and well-mixed planetary boundary layer with elevated lifting condensation level (LCL), forming a negative feedback mechanism (Duveiller et al., 2021; Leung et al., 2024; King et al., 2024). This acceleration is consistent with the biophysical mechanism whereby albedo reduction from idealized afforestation increases sensible heat flux and Bowen ratio, driving stronger buoyant turbulence and the development of a deeper, well-mixed planetary boundary layer (Bright et al., 2017). The resulting enhanced entrainment of dry air and faster rise of the lifting condensation level relative to planetary boundary layer top growth suppress low-level cloud formation, particularly over water-limited mid-to-high latitude land regions (Ek and Holtslag, 2004; Bonan, 2008). In stark contrast, the REAL scenario suppressed the global CLDLOW loss rate (from −0.0153 ± 0.0030 % yr⁻¹ in CTL to −0.0148 ± 0.0029 % yr⁻¹) (Fig. 1e and h). Meanwhile, the BE50 scenario induced a marginal acceleration in the decline (−0.0157 ± 0.0032 % yr⁻¹), equating to a factor of 1.03 globally and 1.23 specifically over land (reaching −0.0207 ± 0.0024 % yr⁻¹) (Fig. S11). Furthermore, all simulations consistently suppressed the loss rate of low-level clouds over the ocean, modifying the trends to −0.0142 ± 0.0039 % yr⁻¹ for REAL and −0.0137 ± 0.0042 % yr⁻¹ for BE50 (Fig. S11). These findings highlight a complex spectrum of biophysical feedback, while extreme land-cover transitions, such as the idealized conversion of grassland or shrubland to forest, trigger strong feedback that accelerate low-level cloud loss, realistic afforestation act as crucial stabilizers against this cloud cover loss.

Figure 1 Simulated long-term trends and seasonal variations of area-weighted global mean cloud cover. Panels display the time series (2015–2100) for total (a–d), low (e–h), medium (i–l), and high-level (m–p) cloud cover under CTL, AF50, BE50, and REAL scenarios. In each panel, heatmaps show the monthly variations, while blue and solid lines track the global annual means. Linear trend lines (dashed) and multi-year means (horizontal solid lines) are superimposed, with statistical significance (p) and regression coefficients (R) provided. The shaded areas indicate the 95 % confidence intervals (CI). Uncertainty propagation for baseline subtractions was performed via quadrature addition of standard errors (Taylor, 1997; Mastrandrea et al., 2010).

Medium-level clouds (CLDMED) exhibited a consistent downward response to land-use forcings, with trends ranging from −0.0101 ± 0.0017 % yr⁻¹ in CTL to −0.0121 ± 0.0019 % yr⁻¹ in AF50, idealized scenarios (AF50 and BE50) induced slightly steeper trends (Fig. 1). Conversely, high-level clouds (CLDHGH) displayed robust positive trends across all land use forcing simulations, with trends increasing from 0.0053 ± 0.0028 % yr⁻¹ in CTL to 0.0082 ± 0.0028 % yr⁻¹ in AF50. The CLDHGH response is associated with large-scale thermodynamic variability; however, the linear models explain only a moderate fraction of the overall variability (Fig. 1m–p). Therefore, these diagnostics do not support assigning the high-cloud response to a single dominant mechanism. Crucially, both idealized and realistic afforestation amplified this upward trend in CLDHGH. Relative to CTL run, AF50 accelerated the global upward trend by a factor of 1.5, an effect that was heavily magnified over land (by a factor of 2.28 ) compared to the ocean (1.41). Similarly, REAL enhanced it by about 1.3 times (2.17 times on land and 1.16 times on ocean). Ocean changes were generally insignificant, possibly due to the buffering effect of remote atmospheric transport (Portmann et al., 2022).

The robust amplification of the CLDHGH positive trend by afforestation, coupled with the acceleration of CLDLOW loss, highlights a critical trade-off in land-based climate mitigation (Luyssaert et al., 2018). While greenhouse gas-induced warming is known to elevate high-cloud tops and amplify climate warming (Norris et al., 2016; Zelinka et al., 2016; Luo et al., 2024a). Our findings reveal a compounding risk, we demonstrate that future large-scale afforestation, particularly under idealized scenarios, could inadvertently exacerbate this warming by reinforcing the positive high-level cloud feedback. Moreover, this effect is further compounded by the pronounced accelerated depletion of low-level clouds over land (as seen in the AF50 scenario), which triggers adverse biophysical feedback that fundamentally compromises the anticipated cooling benefits of forest expansion. Conversely, the REAL scenario exerts a minimal climate perturbation, effectively buffering the background cloud loss. This resilience is largely driven by the spatially targeted allocation of forest expansion (King et al., 2024; Weber et al., 2024; Fig. S12), which constrains cloud feedback within the bounds of natural variability. However, even realistic afforestation carries the risk of exacerbating the increasing trend of high-level clouds, thereby potentially inducing positive radiative forcing at the top of the atmosphere. This indicates that large-scale LUCC affects the lower and upper atmosphere in different ways (Dror and Feingold, 2026). The divergent responses of low- and high-level clouds indicate that their controlling processes may differ. For CLDHGH, surface-flux perturbations may influence deep convection and vertical moisture transport, whereas upper-tropospheric temperature, humidity, and circulation adjustments can also affect the high-cloud response. Because CLDHGH aggregates different high-cloud types and our diagnostics do not separate cirrus from convective anvils, we cannot determine the relative importance of surface-forced convective influences and large-scale dynamical or thermodynamic adjustment. We therefore interpret the CLDHGH increase as an integrated high-cloud response rather than evidence for a uniquely identified dominant pathway (Zhang et al., 2016; King et al., 2024).

3.2 Spatial response of cloud cover to LUCC forcings

Distinct land use forcings exerted pronounced spatially asymmetric impacts on long-term cloud cover trends, characterized by stark geographic heterogeneity between hemispheres and across land-ocean partitions (Figs. 2 and S12–S15). In the Northern Hemisphere (NH), both afforestation and bioenergy expansion triggered widespread negative CLDTOT trends. Notably, BE50 scenario induced the steepest overall decline (mean CLDTOT trend: −0.0045 % yr⁻¹), amplifying the loss rate by a factor of 2.4 compared to the AF50 scenario (mean: −0.0019 % yr⁻¹). However, the most striking spatial anomaly emerged in the CLDLOW regime under AF50 (mean NH trend: −0.0060 % yr⁻¹). Negative CLDLOW trends were strongly localized over mid-to-high latitude landmasses in North America and Eurasia, with the intensity of loss in these land regions being 20 times that of the overall Southern Hemisphere (SH). This phenomenon, indicative of exacerbated seasonal losses during summer, was also evident in BE50 and REAL simulations, suggesting a general latitude dependency (Figs. 3, S12 and S13). Furthermore, CLDMED trends displayed regional divergence, with BE50 inducing CLDMED increases in Australia and Southeast Asia while AF50 amplified losses in those same regions.

Figure 2 Spatial patterns of long-term total and low-level cloud cover trends driven by diverse land-use scenarios. Maps display the anomalous annual trends (% yr⁻¹) of total cloud cover (a–f) and low-level cloud cover (g–l) for both the Northern and Southern Hemispheres. The visualized trends represent the differences (denoted by symbol Δ) between the specific land-use forcings (idealized afforestation, bioenergy expansion, and realistic afforestation) and the control simulation. Blue shading indicates a decreasing trend in cloud cover, while red shading depicts an increasing trend. The cross-hatching highlights grid cells where the long-term trends are statistically significant (p<0.05). Panel headers provide the regional mean trend, alongside the averages of all positive (Pos) and negative (Neg) grid-level trends within the respective domains.

Figure 3 Seasonal and latitudinal distributions of anomalous cloud cover trends under varying land-use scenarios. Heatmaps display the zonally averaged trend differences (% yr⁻¹) across the seasonal cycle (x-axis, months) and latitudes (y-axis, 90° S to 90° N). Columns differentiate the responses for total (a, e, and i), low-level (b, f, and j), medium-level (c, g, and k), and high-level (d, h, and l) cloud fractions. Rows correspond to the idealized afforestation (a–d), bioenergy expansion (e–h), and realistic afforestation (i–l) scenarios relative to the control run. Black asterisks mark grids where the trends are statistically significant at the 95 % confidence level (p<0.05).

Figure 4 Spatiotemporal change of temperature at reference height (TREFHT) and shortwave cloud radiative effect (swCRE). (a–d) Spatial patterns and zonal mean distribution of TREFHT induced by different land use and land cover change. (e–l) Seasonal trends in TREFHT and swCRE grouped by month in different simulation experiments. Given that the global climatological mean of swCRE is approximately −47 W m⁻² (Loeb et al., 2024), an enhanced positive trend in swCRE indicates that shortwave radiative cooling of clouds is either weakening or amplifying warming. The black hatching and stippling denote the grid cells where trends are statistically (p<0.05) and the shade in (d) indicates standard deviation.

At the global scale, the investigated LUCC forcings drove a spatial contraction of the cloud loss footprint, paradoxically accompanied by an intensification of local loss magnitudes (Fig. S16). For instance, while the AF50 scenario contracted the global CLDTOT loss area by 6.01 %, the severity of cloud depletion within these confined regions escalated significantly (Fig. S16a and b). This non-linear response demonstrates that idealized afforestation does not trigger a uniform global reduction in cloudiness; instead, it funnels biophysical perturbations into highly sensitive terrestrial regions, thereby generating localized, intense warming hotspots (Fig. 4a). The formation of such regional warming hotspots is largely associated with areas where grasslands or shrublands have transitioned into forests (Fig. S8). Furthermore, both REAL and BE50 scenarios amplified CLDHGH trend increases over Australia, and REAL induced a more intense mean CLDHGH trend in the SH compared to AF50 (Fig. S12g and h). This implies that even realistic afforestation could inadvertently exacerbate regional warming by accelerating high-cloud formation over specific SH oceanic sectors. Additionally, it is noteworthy that realistic afforestation significantly exacerbates the rate of low-level cloud loss on the western side of South America, even though the forest expansion in this scenario primarily occurs on the southeastern side of the region. This pronounced spatial decoupling highlights critical remote biophysical teleconnections, suggesting that LUCC-induced, cloud-mediated drought risks could be particularly severe in these dynamically linked regions (Li et al., 2025; Kan et al., 2025).

3.3 The response of temperature and cloud radiative effect to LUCC

Beyond the direct (local) surface albedo effect, distinct land use forcings influence temperature trends and their spatial distribution through indirect (non-local) cloud radiative effects (Schneider, 1972). While increases in temperature at reference height (TREFHT) are ubiquitous across all simulations (Fig. 4e–h), idealized afforestation does not merely amplify this background warming; rather, it coalesces these thermal trends into highly concentrated regional hotspots (Figs. 4a and S17a, e). This localized warming is predominantly driven by the mid- to high latitudes of the Northern Hemisphere. To quantify this spatial aggregation, we applied the Gini coefficient to the positive temperature trends in the Northern Hemisphere. Notably, the coefficient shifts from 0.25 in the CTL simulation to 0.31 under the AF50 scenario (Fig. S17a and e). This stark exacerbation in spatial inequality demonstrates that large-scale idealized afforestation induces a highly skewed, uneven warming pattern – disproportionately penalizing specific latitudinal bands rather than distributing the thermal burden uniformly. The pronounced spatial aggregation of warming induced by the AF50 scenario is fundamentally driven by the extensive conversion of native grasslands and shrublands to forests (Fig. S8c). Because these non-forest biomes are inherently concentrated across mid-to-high latitudes (Olson et al., 2001; Dixon et al., 2014), their large-scale replacement disproportionately localizes the adverse biophysical warming effects within these specific latitudinal bands. The expansion of bioenergy crops substantially redistributes the warming hotspots induced by AF50, the intensity is alleviated in many regions (e.g., Southern Ocean, Antarctica), while opposite changes occur elsewhere, such as the Western Pacific, where BE50 reverses AF50's cooling to a warming trend. A similar redistribution of hotspots is observed when comparing AF50 with realistic afforestation; notably, REAL mitigates or offsets the significant warming trends seen in AF50 across mid-to-high latitudes of North America and Eurasia.

Figure 5 Long-term trends of key physical terms (decomposed by the surface energy balance equation) induced by different LUCC. T_a represents the atmospheric feedback term, albedo denotes the surface albedo term, and SWdown (SW_↓) and LWdown (LW_↓) indicate the downward shortwave and longwave radiation terms, respectively. The black hatch indicates the grids that passed the test of significance (p<0.05).

The zonal mean responses of TREFHT exhibit stark divergences (Fig. 4d). In the mid-to-high latitudes of the NH, the idealized afforestation drives substantial warming, in sharp contrast to the pronounced cooling induced by the REAL scenario. This boreal warming is intricately coupled with severe cloud cover loss; specifically, the loss of low-level clouds triggers an anomalous positive swCRE trends that substantially amplifies surface heating (Figs. 4, 5, and S18–S21). Crucially, the top-of-atmosphere (TOA) shortwave cloud radiative effect (swCRE) serves as a robust diagnostic, vividly capturing this signature of cloud-mediated warming amplification (Figs. S20 and 21). The regions with larger cloud loss magnitudes and stronger temperature rise trends largely overlap, particularly in the boreal region (50–90° N). Relative to the global average, AF50 and BE50 significantly amplified warming here (factors of 2.15 times and 1.71 times, respectively), while REAL induced cooling (factor of −1.51 times). This pronounced warming amplification phenomenon is primarily driven by an enhanced atmospheric feedback trend, further amplified by the accelerated reduction of clouds in the middle and low-level clouds (Fig. 5). AF50 caused the low-level cloud loss rate in this region to be 8.80 times the global average, and BE50's relative amplification was even higher (13.63 times). Furthermore, the total cloud cover decline significantly contributed to the temperature rise, showing 8.88 times the global level in AF50 and by a factor of 26.12 in BE50. Therefore, under the idealized LUCC scenario, the significant warming phenomenon simulated in the boreal region is partly attributed to the enhanced positive trend of swCRE caused by mid-to-low-level cloud loss, although atmospheric feedback plays a dominant role. A significantly enhanced trend in swCRE is simulated on the western side of South America, which is associated with accelerated loss of low-level clouds. Furthermore, at the global scale, the sensitivity of swCRE to temperature is attenuated under idealized afforestation and bioenergy expansion scenarios relative to the control. Specifically, the response weakens from −0.746 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$ in the control to −0.741 and −0.742 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$ in the AF50 and BE50 experiments, respectively. Conversely, realistic afforestation exhibits an opposing trend (−0.749 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$). In other words, the negative shortwave cloud feedback is weakened by approximately 0.67 % and 0.54 % under the two idealized scenarios, whereas realistic forest expansion strengthens this negative feedback by ∼ 0.40 % (Fig. S22).

3.4 Forest and bioenergy-induced hemispheric cloud trends asymmetry

We used the Gini index to quantitatively assess the inequality and spatial concentration of global cloud cover trends due to LUCC (Zhao et al., 2025). A Gini index closer to one signifies a highly unequal distribution, with trends strongly concentrated in specific latitude bands. For CLDTOT, relative to the CTL, idealized afforestation significantly intensified concentrated cloud losses in both hemispheres (e.g., SH downward trends Gini index rose from 0.48 to 0.52; Figs. 6 and S23d), suggesting amplified local land use impacts (Hua et al., 2023). This concentration effect was most dramatic for CLDLOW in the NH, where AF50 caused the Gini index to rise sharply from 0.41 (CTL) to 0.50 (Figs. 6f and S23f), indicating a greater concentration of cloud decline areas. Bioenergy expansion induced milder asymmetry overall, but still concentrated SH increases in specific CLDTOT and CLDHGH hotspots (Fig. S24).

Figure 6 Asymmetry of cloud trends between hemispheres. Lorenz curves of long-term CLDTOT (a–d), CLDLOW (e–h), CLDMED (i–l), CLDHGH (m–p) trend in the Northern (0 to 90° N) and Southern Hemispheres (0 to 90° S) during the 2015–2100 in CTL run. The share indicates the percentage of an increase or decrease in the trend. The results of the AF50, BE50, and REAL experiments can be found in Figs. S23–S25.

In stark contrast to idealized expansion, realistic afforestation consistently exerted the minimal influence on cloud trend heterogeneity, remaining closest to the CTL simulation (Fig. S25). For CLDTOT, REAL scenario slightly broadened NH increasing trends (Gini index from 0.54 to 0.50) with limited alteration to the overall spatial patterns. For CLDLOW, REAL scenario increased the concentration of NH downward trends, but its impact was far less pronounced than that of AF50 (Gini index from 0.41 to 0.43), this means that realistic afforestation will concentrate the loss of low-level clouds in specific areas. Collectively, these findings reveal that idealized LUCC forcings profoundly disrupt inter-hemispheric cloud homogeneity; by triggering highly concentrated low-level cloud depletion, they risk severe localized climatic perturbations, a hazard largely avoided under the more spatially diffuse impacts of realistic scenarios. Previous studies has established that large-scale land-use changes profoundly perturb atmospheric dynamics over a radius extending at least 200 km (Chambers and Artaxo, 2017; Khanna et al., 2017). Crucially, within this expansive dynamic footprint, the intense spatial aggregation of low-level cloud loss driven by the widespread conversion of native grasslands and shrublands into forests, holds the potential to precipitate severe regional droughts. Analogous to the hydrological extremes induced by highly skewed precipitation patterns (Lesk and Mankin, 2026), this geographically concentrated cloud loss disrupts the local energy and moisture balance, acting as a potent catalyst for cloud-mediated drought risks (Miralles et al., 2014; Zhou et al., 2019).

3.5 Machine learning identifies influential predictors of LUCC-induced cloud trends

Having established that distinct LUCC scenarios markedly alter the spatial distribution of global cloud-cover trends, we employed an interpretable machine-learning framework to identify the most influential predictors of these changes and characterize their relative importance and model-learned dependencies.

Figure 7 Ranked drivers to cloud cover trends induced by different land-use change scenarios, as revealed by the machine learning model. The SHAP (SHapley Additive exPlanations) summary plots illustrating the impact of each climate features on the predicted changes in (a–c) total cloud (CLDTOT), (d–f) low-level cloud (CLDLOW), (g–i) medium-level cloud (CLDMED), and (j–l) high-level cloud (CLDHGH) under the AF50, BE50, and REAL scenarios relative to the CTL run. Δ present the difference between sensitivity experiments and CTL. In other words, this represents a machine-learning-driven mechanistic modeling of LUCC-induced cloud cover and climate trends. Each dot represents an individual sample. The x-axis (SHAP value) quantifies the directional impact of a specific feature on the model's predicted cloud trend (i.e., a positive SHAP value drives an increasing cloud cover trend, while a negative value drives a decreasing trend). The right-side color bar signifies the actual physical magnitude of the climate feature itself, transitioning from low values (orange) to high values (purple). The relationship between the color and the x-axis reveals the physical driving mechanism. For instance, if a cluster of purple dots (high feature value) is distributed on the positive side of the x-axis (SHAP > 0), it physically indicates that an increase in this specific climate feature induced by LUCC strongly drives an increasing trend in cloud cover. Conversely, orange dots on the negative side indicate that a decrease in the feature value contributes to a decreasing cloud cover trend. PRECT, total precipitation ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$); TMQ, total precipitable water ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$); TREFHT, air temperature at reference height (K yr⁻¹); TREFHTMX, maximum temperature at reference height (K yr⁻¹); RHREFHT, relative humidity at reference height (% yr⁻¹); QFLX, evapotranspiration ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$); CAPE, convective available potential energy ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$); PBLH, planetary boundary layer height (m yr⁻¹); TREFHTMN, minimum temperature at reference height (K yr⁻¹); SHFLX, sensible heat flux ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$); LHFLX, latent heat flux ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$).

Figure 8 Normalized feature importance for cloud-cover trends under different land-use change scenarios. Δ present the difference between sensitivity experiments and CTL. Feature importance is calculated from SHAP values and normalized across all predictors within each machine-learning model. The resulting percentages quantify the relative importance of predictors to model predictions (Berdugo et al., 2022).

Figure 9 Dependence plots of the top five important variables for total cloud cover (CLDTOT), as revealed by the machine learning model. The first to third rows represent AF50, BE50 and REAL experiments, respectively. The background scatter displays the full samples colored by feature values. The SHAP dependence trend (black solid line) was visualized using Locally Weighted Scatterplot Smoothing (LOWESS) (Cleveland, 1979) with a smoothing fraction of frac = 0.3 (using 30 % of nearest neighbors for local regression) and it = 3 robust iterations to reduce outlier influence, and gray shaded band indicates 95 % confidence interval (CI) of LOWESS regression (calculated by a random subsample of 10 000 grid cells). The corresponding results for CLDLOW, CLDMED, and CLDHGH are presented in Figs. S26, S28, and S30, respectively.

Figure 10 Spatial drivers and contribution of CLDTOT revealed by machine learning response to different land use change. The first to third rows represent AF50, BE50 and REAL experiments, respectively. The contribution is expressed in SHapley Additive exPlanations (SHAP) value (Berdugo et al., 2022). The corresponding results for CLDLOW, CLDMED, and CLDHGH are presented in Figs. S27, S29, and S31, respectively.

Figure 11 Vertical profiles of temperature and relative humidity trends induced by different LUCC. The black asterisks denote grid cells with statistically significant trends (p<0.05).

Across all simulations, the machine-learning analysis identifies changes in total precipitable water (ΔTMQ) as the highest-ranked predictor of CLDTOT trends, with normalized SHAP importance values of 20.5 %, 26.3 %, and 23.1 % under the AF50, BE50, and REAL scenarios, respectively (Figs. 7a–c and 8a–c). The SHAP dependence relationships between TMQ and CLDTOT are consistently positive (Fig. 9a, f, and k), indicating that lower TMQ values are associated with more negative model-predicted CLDTOT trends, whereas higher TMQ values are associated with more positive predictions, which indicates that large-scale LUCC disturbances to the hydrological cycle are the primary macroscopic factor controlling long-term changes in total cloud cover (Pielke et al., 2007). Notably, the dependency of CLDTOT on TMQ under idealized afforestation is highly linear and significantly stronger than in other scenarios, the particularly strong and approximately linear SHAP dependence under AF50 indicates that TMQ is especially influential in the model predictions for this scenario (Fig. 9a, f, and k). Spatially, however, the regulatory role of TMQ varies by location; the TMQ increase trend induced by idealized forest expansion mainly contributes to the CLDTOT increase trend in tropical regions, while outside the tropics, it is primarily negative contributions (Fig. 10a). These results indicate that the rapid CLDTOT decline across Northern Hemisphere land (Fig. 2) is fundamentally governed by a chronic deficit in atmospheric moisture supply, and this manifests specifically as a trend of significant loss in relative humidity (Fig. S32b). These dynamics corroborate previous modeling evidence demonstrating that extensive midlatitude land-cover transitions (from C₃ grasslands to forests) severely exacerbate lower-tropospheric moisture deficits, driving a concomitant and pronounced depletion of cloud cover (Laguë and Swann, 2016). Particularly in water-limited regions of the mid- to high-latitudes of the Northern Hemisphere, the transpiration induced by forests replacing grasslands tends to be weak, which implies that the efficiency of water vapor transport from forests to the atmosphere is lower (Breil et al., 2021). The temperature rise and relative humidity decrease caused by idealized forest expansion are not limited to the surface. Simulations show that the warming and drying trends extend vertically through the troposphere (Fig. 11). The decrease in relative humidity causes the lifting condensation level (LCL) of water vapor to rise. If the elevation of the LCL exceeds the deepening of the planetary boundary layer, boundary-layer moisture fails to reach saturation, thereby inhibiting cloud development or forcing emerging clouds to prematurely dissipate in the overlying subsaturated air (Romps, 2017; Kutta and Hubbart, 2023). This thermodynamic decoupling elucidates how idealized forest expansion drives a synchronized depletion of mid-level clouds (accelerating by a factor of 1.25 over land). Specifically, surface thermodynamic stress dampens vertical moisture transport by altering boundary-layer stability and depth, effectively suppressing mid-level cloud development. Conversely, the REAL scenario averts this desiccation by optimizing surface-atmosphere heat partitioning; it sustains robust evapotranspiration efficiency and promotes atmospheric cooling, thereby effectively neutralizing the albedo-driven warming (Figs. 11 and S34). Previous modeling studies have shown that progressive afforestation in mid-latitudes may suppress moisture supply, thereby causing tropospheric drying; the reduction in cloud fraction across vertical levels is precisely the result of this atmospheric response (Laguë and Swann, 2016). Therefore, in the idealized afforestation simulation, relative humidity is the second-ranked predictor of CLDTOT trends, with a normalized SHAP importance of 14.1 % (Fig. 8a). This total precipitable water-dominated attribution broadly extends to regional total cloud cover losses under bioenergy expansion and realistic afforestation. However, the subsequent hierarchy of contributing factors exhibits distinct scenario-dependent variations, reflecting varying interplay between relative humidity and surface temperature extremes (Fig. 8b and c).

For CLDLOW, relative humidity is the highest-ranked predictor under all three LUCC scenarios, with normalized SHAP importance values of 20.4 %, 20.4 %, and 15.3 % for AF50, BE50, and REAL, respectively (Figs. 7d–f and  8d–f). Under AF50, decreasing relative humidity is associated with increasingly negative SHAP values for CLDLOW, with the strongest spatial correspondence concentrated over the mid-to-high latitudes of the Northern Hemisphere (Fig. S27a). In these regions, idealized afforestation causes surface darkening and a substantial reduction in albedo, leading to increased absorption of surface net solar radiation (Figs. S35–S37), particularly during the late afforestation stage (2081–2100), idealized forest expansion induces a pronounced decrease in land surface albedo by 0.0142, a signal that is approximately 2.4 times larger than the global mean (−0.0058). This albedo anomaly triggers a concomitant increase in net surface radiation over land by 2.43 W m⁻², exceeding the global average (0.74 W m⁻²) by a factor of roughly 3.3 (Fig. S37). The additional energy is preferentially partitioned into sensible heat flux (elevated Bowen ratio), suppressing the increase in latent heat flux. Over land areas, the sensible heat flux increases by 0.98 W m⁻², whereas the latent heat flux increases by only 0.92 W m⁻², generating strong buoyant turbulence that promotes the development of a deeper and well-mixed planetary boundary layer. Furthermore, this implies that as larger-scale grasslands or shrublands are converted into forests, while the initial latent-heat-dominated evaporative enhancement may slightly promote cloud cover, the subsequent substantial escalation in sensible heat flux significantly accelerates the loss of low- and mid-level clouds. Enhanced vertical mixing and entrainment of drier free-tropospheric air dilute boundary-layer moisture while warming the lower troposphere, jointly driving the pronounced decline in relative humidity and the rise in maximum temperature simulated under the idealized afforestation scenario (Fig. S11; Spracklen and Garcia-Carreras, 2015). At the same time, idealized afforestation induces a significant increase in temperature in response to the substantial and rapid loss of CLDLOW (Fig. S32d, e, and g), which indicates that the relative humidity decline, and temperature rise trends induced by idealized afforestation jointly amplify the loss rate of CLDLOW, and it exhibits a strong linear dependency (Fig. S26a and b). As established by prior studies, afforestation across the NH mid-to-high latitudes induces profound surface darkening, which sharply reduces albedo and amplifies solar radiation absorption, ultimately driving surface warming (Meissner et al., 2003; Swann et al., 2010). Crucially, the magnitude of this albedo decline and the subsequent positive radiative forcing per unit of afforested area are disproportionately greater in these boreal domains than in any other region (Foley et al., 1994; Weber et al., 2024). The simulation proves that the regions with rapid decline trends in CLDLOW induced by idealized afforestation simultaneously experience rapid rise trends in temperature and rapid decline in relative humidity (Figs. 2 and S32b, d, e, and g). While albedo reduction acts as the initial biophysical driver of local surface warming across NH land areas (Arora and Montenegro, 2011), it cannot directly precipitate the long-term decline of cloud cover. As revealed by our machine learning attribution, the critical missing bridge is governed by atmospheric moisture dynamics. Specifically, under idealized afforestation, the albedo-induced surface warming fundamentally alters surface energy partitioning, heavily skewing it toward sensible heat fluxes. This intense sensible heating drives the development of a deeper, vigorously mixed planetary boundary layer (Fig. S32f). Consequently, enhanced dry-air entrainment from above, coupled with a lifting condensation level that ascends faster than the planetary boundary layer top, severely depresses relative humidity (Fig. S32b). This thermodynamic cascade provides the definitive physical link between surface albedo reduction and low-level cloud cover loss. This is highly consistent with our machine learning results, which parse relative humidity, closely followed by maximum temperature, as the paramount driving factors for accelerated low-level cloud cover loss.

Machine learning identifies different LUCC-induced trends in total precipitation (ΔPRECT) and atmospheric total precipitable water (ΔTMQ) as the paramount drivers governing mid- and high-level cloud dynamics (Fig. 7g–l), in particular, these two factors contribute more than 40 % to the changing trends of mid-level clouds (Fig. 8g–i). Large-scale idealized afforestation causes insufficient water supply in the mid-to-high latitude regions of the NH, which also inhibits the development of CLDMED, but the impacts of bioenergy expansion and realistic afforestation are overall positive (Fig. S29). More critically, high-level clouds can exert a net warming effect on the climate system through their strong longwave greenhouse effect (Hartmann et al., 1992; Stephens, 2005). In the idealized afforestation experiment, the simulated increase in CLDHGH over the mid-to-high latitudes of the Northern Hemisphere coincides with upper-tropospheric cooling and increased relative humidity (Fig. 11a and d), indicating that the high-cloud response is closely associated with changes in the upper-tropospheric thermodynamic environment. However, because CLDHGH is an aggregate pressure-based diagnostic, these relationships do not establish whether the simulated increase primarily reflects cirrus, convective anvils, surface-flux influences on deep convection, or indirect large-scale dynamical adjustment. The concurrent increase in the top-of-atmosphere longwave cloud radiative effect (Fig. S20b) demonstrates the radiative consequence of the enhanced high-cloud cover, but does not distinguish among these formation pathways. Specifically, in the Northern Hemisphere, the mean positive trend of TOA swCRE rises from 0.0465 W m⁻² in the CTL run to 0.0490 W m⁻² in AF50, reaching 0.0511 W m⁻² in BE50 and 0.0531 W m⁻² in the realistic afforestation scenario. Similarly, regarding the TOA longwave cloud radiative effect in the Northern Hemisphere, the mean positive value surges from 0.0099 W m⁻² in the CTL run to 0.0155 W m⁻² in AF50. Therefore, idealized afforestation not only accelerates warming by accelerating the loss of mid- and low-level clouds, but also amplifies warming by promoting the increase of high-level clouds, and the warming induced by high-level clouds also exists in the bioenergy expansion and realistic afforestation simulations. However, in contrast, afforestation in real-world scenarios effectively optimizes the distribution of heat between the surface and the atmosphere, allowing atmospheric cooling and evapotranspiration efficiency (a positive manifestation of hydrothermal coupling) to dominate, partially offsetting the warming effects caused by decreased albedo. Admittedly, under realistic afforestation scenarios, the pronounced cooling simulated across the Northern Hemisphere mid-to-high latitudes may not be driven exclusively by positive cloud cover anomalies, as non-local temperature advection likely also plays a substantial modulating role (Kan et al., 2025).

4 Discussions

Our simulation results indicate that two distinctly different future land use methods may significantly alter long-term cloud cover trends in terms of their impact on the magnitude of climate warming, particularly large-scale idealized forest expansion exacerbates cloud-mediated climate warming by substantially accelerating the loss of mid-low-level clouds and the increase of high-level clouds. Bioenergy expansion also imposes negative impacts on cloud cover, but the overall intensity is weaker than that of idealized afforestation; this relatively mild negative impact may stem from lower LAI and reduced transpiration, as in previous assessments of this scheme, where bioenergy expansion may be accompanied by depletion of soil moisture and hinder the accumulation of atmospheric humidity (Le et al., 2011; Cheng et al., 2024). We believe that the negative impacts of these two large-scale idealized LUCC forcings on cloud cover originate from their own expansion areas, because in our simulations, the land for forest expansion is taken from grasslands or shrubs, which may lead to an increase in surface radiation absorption; the increased radiation warms the bottom of the convection layer in the form of sensible heat, thereby imposing more stringent conditions for cloud formation. In contrast, in the bioenergy expansion simulation (irrigated and rained sugarcane), the land comes from the land of other bioenergy crops, which means a reduction in the area of other crops with higher albedo, for example, switchgrass has an overall higher albedo than sugarcane and thus has greater climate cooling potential (Abraha et al., 2021; Lei et al., 2023). When other possible bioenergy crops with greater climate cooling potential are replaced by sugarcane, it may suppress cloud formation due to insufficient cooling, thereby inducing mild cloud loss. Furthermore, since our idealized afforestation scenarios involve converting grasslands or shrubs to forests, although this conforms to the Food and Agriculture Organization (FAO) definition (FAO, 2012), meaning our simulations are neither forest enhancement nor reforestation, it is also essential to consider the potential ecological consequences of such land type conversions. Grasslands, as one of the important terrestrial ecosystems, support unique biodiversity but are also relatively fragile; afforestation in these areas may lead to the loss of flora and fauna, reduced water availability, and increased fires (Veldman et al., 2015a, b; Parr et al., 2024). Therefore, although the primary goal of idealized afforestation is to enhance terrestrial carbon absorption, it carries the risk of ecological degradation (Cao, 2008; Xiao et al., 2020). Here, our idealized forest expansion simulation experiment is not intended to suggest implementing such forest expansion actions; on the contrary, the signal we are sending is cautionary.

On the other hand, we found that idealized and realistic afforestation suppressed the loss trends of total and low-level cloud over the oceans, but idealized bioenergy expansion accelerated their loss trends (Fig. S11). The question we need to answer here is actually through what pathways LUCC disturbances on land propagate energy to the oceans and influence cloud trends. According to previous simulation studies on deforestation, changes in albedo are considered the fundamental source of ocean responses; deforestation increases surface albedo, thereby reducing energy absorption at the surface, so no more heat is allocated to sensible heat flux to heat the convection layer, this process leads to cooling of the convection layer over the oceans (Davin and de Noblet-Ducoudré, 2010). In contrast, in our simulations, the albedo reduction caused by forest expansion (Figs. S35–S37) induces the surface to absorb more solar radiation; the extra radiation is redistributed in the form of sensible heat flux to release into the atmosphere and heat the atmosphere, and the increased heat warms the ocean surface with ample moisture supply, making more water vapor condense into clouds. Therefore, forest expansion overall suppresses the loss rate of total cloud and low-level cloud cover in ocean areas (Liu et al., 2023).

Our hypothesis can be verified from Figs. S32 and S33, namely that forest expansion induces increasing trends in maximum temperature and total atmospheric precipitable water over the oceans. However, in bioenergy expansion, we observe overall opposite results (Fig. S32), namely that bioenergy, due to relatively high albedo, does not lead to substantial warming of the oceans, thus there is not enough heat to make moisture evaporate and condense into clouds. In contrast, realistic afforestation overall does not render clouds detrimental to the climate; this difference largely stems from the type of land conversion and the regions involved. In idealized afforestation, forests are derived from grasslands or shrubs, with large amounts of C₃ grasslands distributed in mid-to-high latitude regions of the Northern Hemisphere and C₄ grasslands in tropical areas, whereas forest expansion in realistic afforestation simulations is not like this (Wei et al., 2014; Fig. S8).

Although our simulations show that large-scale idealized afforestation and bioenergy expansion induce significant declining trends in mid- and low-level cloud cover, the SHAP-based spatial attribution from machine learning clearly demonstrates that these two distinct land-use changes actually trigger multiple competing land–atmosphere coupling processes. The net cloud response is essentially the outcome of the interplay among these competing drivers (Figs. S38–S40), even though the dominant factors identified by machine learning differ across scenarios (Fig. 7). The spatial competing patterns, visually represented by the stippled regions, arise from the competition between surface fluxes and boundary-layer dynamics. On one hand, forest expansion increases leaf area index and surface roughness, thereby enhancing turbulent mixing and evapotranspiration (latent heat flux). This process tends to moisten the lower atmosphere and favor cloud formation, as supported by the increasing trend in total precipitable water induced by idealized afforestation (Fig. S32a). On the other hand, the large-scale replacement of grasslands or shrublands by forests inevitably causes a substantial reduction in surface albedo. The resulting increase in absorbed solar radiation is directly manifested as enhanced sensible heat flux (Fig. S32h). This strong sensible heating deepens the planetary boundary layer and promotes the entrainment of dry air from aloft (as indicated by the rising trend in boundary layer height), generating a strong cloud-suppressing effect. Specifically, in water-limited regions of the Northern Hemisphere mid-to-high latitudes, the modest increase in transpiration resulting from the conversion of grasslands/shrublands to forests is insufficient to offset the strong drying effect driven by sensible heat flux. Compared with the substantial decline in relative humidity, the increase in total precipitable water remains relatively small. This leads to intense spatial competition between near-surface relative humidity (RHREFHT) and maximum temperature at reference height (TREFHTMX), ultimately causing a sharp reduction in low-level clouds (Fig. S38d and e). In the subtropical regions under the BE50 scenario, due to the relatively low transpiration efficiency of bioenergy crops, the dynamic drying effect associated with the rise in planetary boundary layer height (PBLH) becomes the more dominant factor (Fig. S39e). In contrast, the more realistic REAL scenario is constrained by the discontinuous spatial distribution of afforestation patches. Local advective exchange within microclimates exerts a strong physical buffering effect, resulting in a markedly fragmented pattern of large-scale competition. This fragmentation helps mitigate the widespread cloud loss observed under the idealized scenarios (Fig. S40d–f).

These findings align with recent studies highlighting the regime-dependent nature of cloud responses to land-cover change. For example, Xu et al. (2022) and Leung et al. (2024) demonstrated that the sign and magnitude of forest effects on clouds vary substantially with background climate, soil moisture availability, and aerosol conditions. Our results reinforce the view that large-scale idealized afforestation does not produce universally positive cloud feedback; instead, the net cloud response is determined by the balance between competing local processes, which is particularly unfavorable for low-cloud maintenance in water-limited mid-to-high latitude regions.

The high-level cloud response requires a different interpretation from the boundary-layer mechanism identified for low clouds. Because CLDHGH is defined by cloud-top pressure above 400 hPa, it represents an aggregate high-cloud diagnostic rather than a specific cloud type. Changes in surface albedo and turbulent heat fluxes can alter boundary-layer buoyancy and may influence deep convection and vertical moisture transport, providing a potential pathway from land-surface forcing to the upper troposphere. At the same time, the simulated upper-tropospheric cooling and moistening (Fig. 11), together with the diagnosed large-scale atmospheric adjustment, indicate that the thermodynamic and dynamical environment in which high clouds form is also modified by LUCC. These pathways are not mutually exclusive: surface forcing can act as an upstream perturbation, while convective transport and large-scale atmospheric adjustment can mediate the subsequent upper-tropospheric response. However, because the present diagnostics do not distinguish cirrus from convective anvils or quantitatively isolate the contributions of these pathways, we do not attribute the CLDHGH increase to any single dominant mechanism. Instead, we interpret it as an integrated response to changes in the upper-tropospheric thermodynamic and dynamical environment.

Recent observational studies have increasingly highlighted that cloud responses can fundamentally alter the net biophysical climate forcing of land-use changes, though the direction of this feedback is highly regime-dependent. For example, recent quantifications of forest loss present contrasting cloud radiative adjustments: Luo et al. (2024b) demonstrated that global deforestation typically reduces cloud cover, thereby introducing a warming effect that partially offsets the cooling generated by increased surface albedo. Conversely, Dror and Feingold. (2026) observed that deforestation in the Amazon basin increases low-level cloudiness, which acts synergistically to roughly double the top-of-atmosphere cooling effect of surface brightening alone. Despite these regional divergences in deforestation impacts, both studies underscore a critical biophysical framework: the interplay between surface albedo and cloud albedo adjustments. Viewed through this framework, our findings for large-scale idealized afforestation reveal a compounding climate penalty. In our AF50 scenario, forest expansion (large-scale conversion of grasslands or shrublands into forests) substantially decreases surface albedo (a direct warming effect), while the concurrent loss of low-level clouds further diminishes shortwave cloud cooling. Consequently, rather than counteracting one another, these two biophysical responses act synergistically to strongly amplify regional climate warming hotspots, particularly in the water-limited boreal regions.

Our findings underscore that large-scale afforestation, despite being widely promoted as a primary climate mitigation strategy, does not universally guarantee a net cooling effect. The exact climatic outcome is not only contingent upon the background climate, but is also heavily governed by the specific pathways of land-cover types conversion. Particularly over the past few decades, there has been a proliferation of afforestation initiatives implemented in regions naturally occupied by grasslands or shrublands (Veldman et al., 2015a; Cao et al., 2011). Furthermore, efforts aimed at enhancing canopy density within existing forests are conceptually distinct from true afforestation; conflating these practices can severely distort objective assessments of the true climate mitigation efficiency of forest expansion. It is imperative to reiterate that the idealized afforestation scenario simulated in this study strictly aligns with the standard definition established by the FAO. Consequently, our modeling results offer a critical caveat: when large-scale land-cover conversion adheres rigidly to this foundational definition, the climate mitigation potential traditionally attributed to practical afforestation may be substantially overestimated. This implies that anchoring future climate stabilization strategies heavily upon widespread forest expansion pathways introduces profound risks and uncertainties (Arora and Montenegro, 2011; Veldman et al., 2015a; Weber et al., 2024).

Importantly, the climatic impacts of large-scale LUCC are not strictly confined to the regions where land-cover conversion occurs. Localized biophysical perturbations can instead cascade into far-reaching non-local climate teleconnections (Portmann et al., 2022; Hua et al., 2023). Several methods have been developed to disentangle local and non-local biophysical responses to LUCC (Zhao et al., 2026). In particular, checkerboard experiments can be applied to both land and atmospheric variables without variable-specific physical assumptions, and Chen et al. (2022) demonstrated their application to the separation of low-, middle-, and high-level cloud-fraction responses. Because the transient AF50, BE50, and REAL experiments in the present study were not designed for checkerboard separation and we do not apply another dedicated cloud-separation framework here, we do not quantitatively partition cloud-cover changes into local and non-local components. Therefore, the following temperature-advection and wave-activity-flux diagnostics are interpreted as evidence for remote dynamical modulation and teleconnection pathways rather than as a quantitative local or non-local decomposition of the cloud response. Our diagnostic analysis of horizontal temperature advection in the lower troposphere (850 hPa) demonstrates that, under idealized afforestation and bioenergy expansion scenarios, intense local sensible heating induced by reduced surface albedo generates anomalous wind fields that redistribute thermal energy to adjacent regions (Fig. S41). Specifically, in the idealized afforestation (AF50) scenario, pronounced positive anomalies in total temperature advection develop over the Northern Hemisphere mid- and high latitudes, effectively transporting excess surface heat from afforested areas to downwind non-afforested regions. This advective transport serves as a dynamical bridge, explaining why warming hotspots frequently extend spatially beyond the precise footprint of the land-use forcing. In contrast, realistic afforestation (REAL) mitigates excessive local sensible heat accumulation; the associated weakening of westerlies reduces warm advection toward Eurasia, thereby accounting for the regional cooling and limiting the extent of non-local warming.

These non-local teleconnections are further amplified through planetary-wave dynamics in the mid- to upper troposphere. The strong surface thermal anomalies triggered by large-scale idealized LUCC perturb the overlying atmosphere, acting as a deep heat source that excites quasi-stationary Rossby waves. Evaluation of the two-dimensional Plumb wave activity flux at 500 hPa reveals that planetary-wave energy emanating from these thermal perturbations exhibits pronounced horizontal propagation (Fig. S42). In the AF50 scenario, the wave activity flux patterns illustrate clear eastward and poleward energy propagation, which conveys additional energy into non-afforested regions across the Northern Hemisphere mid- and high latitudes. The wave-activity-flux patterns are consistent with a remote dynamical pathway through which atmospheric energy propagation can modify subsidence, static stability, and humidity in regions beyond the land-cover perturbation, thereby modulating cloud cover. In the AF50 scenario, the eastward and poleward propagation of wave activity coincides spatially with remote warming, atmospheric moisture deficits, and relative-humidity reductions in the Northern Hemisphere mid- and high latitudes. These relationships support the presence of non-local atmospheric adjustment, but the present diagnostics do not determine its quantitative contribution relative to local processes. Therefore, the pronounced mid- and low-level cloud loss under idealized afforestation is interpreted as being influenced by local thermodynamic and cloud-radiative processes together with non-local dynamical adjustment, without assigning dominance to either pathway. The planetary-scale wave response to idealized bioenergy expansion is comparatively weaker but similarly promotes non-local warming and modestly amplifies mid- and low-level cloud loss. In stark contrast, realistic afforestation strengthens northerly flows while weakening westerlies over Eurasia, suppressing regional warming and producing a net remote cooling effect (Figs. S41 and S42). This process improves the thermodynamic conditions required for cloud formation, consistent with the simulated increases in mid- and low-level cloud cover over the region (Fig. 2).

Due to the intrinsic uncertainties inherent in individual climate models, particularly concerning the simulation of cloud feedbacks (Ceppi et al., 2024), the evaluation of long-term impacts arising from external forcings such as land use and land cover change can often be robustly achieved through well-validated single-model simulations, as demonstrated in numerous prior studies (Snyder et al., 2004; Arora and Montenegro, 2011; Bauer et al., 2025). However, considering that CESM may be overly sensitive to aerosol-cloud interactions or thermal disturbances, the simulated cloud cover trend changes induced by LUCC may be an upper limit. While model ensembles or multi-model intercomparison frameworks have been employed to enhance confidence in conclusions (Falloon et al., 2014), they introduce limitations such as computational expense, potential overconfidence due to model dependence – where shared code, parameterizations, and biases among models undermine the assumption of independence (Knutti et al., 2010; Abramowitz et al., 2019), and challenges in combining projections from a small or unevenly distributed set of models (von Trentini et al., 2020). Despite this, we still recommend, where computational resources permit, the integration of models with different equilibrium climate sensitivity (ECS) characteristics to adequately distinguish between signal and model structure biases (Khan et al., 2025).

Cloud parameterization uncertainties in CESM warrant particular scrutiny because they directly modulate both the local biophysical pathways and the non-local dynamical teleconnections through which LUCC influences cloud cover. CESM employs the Community Atmosphere Model version 6 (CAM6), whose stratiform cloud microphysics is governed by the MG2 two-moment scheme (Morrison et al., 2020). This scheme prognoses mass and number concentrations for cloud liquid, ice, rain, and snow and incorporates an aerosol-aware heterogeneous ice nucleation parameterization (Fan et al., 2016). However, several aspects of this framework introduce substantial uncertainty. A documented issue with the cloud ice number limiter (nimax) in the default CESM configuration inadvertently suppressed ice production in mixed-phase clouds, altering global cloud phase partitioning and contributing to the model's high ECS (Gettelman et al., 2019; Zhu et al., 2022). Subsequent analyses have shown that ice nucleation processes influence cloud feedback across a wide range of tropospheric levels and latitudes, not only in extratropical low clouds, and that correcting or perturbing these processes can substantially modify shortwave cloud feedback (McGraw et al., 2023; Zhu et al., 2022). In addition, the representation of sub-grid cloud variability in warm rain formation relies on an enhancement factor for auto-conversion that is derived from the Cloud Layers Unified by Binormals (CLUBB) turbulence scheme. This factor has been shown to be substantially overestimated over much of the ocean, leading to excessive auto conversion rates, overly strong aerosol–cloud interaction sensitivity, and biases in liquid water path and cloud optical thickness (Wang et al., 2022). Perturbed-parameter ensemble experiments further demonstrate that simulated cloud properties, precipitation, and radiative effects in CAM6 are highly sensitive to individual microphysical parameters, including those controlling auto conversion and accretion rates, ice fall speed, and the ice-to-snow conversion threshold. These sensitivities are especially relevant for LUCC experiments because afforestation and bioenergy expansion modify surface sensible and latent heat fluxes, planetary boundary layer stability, and near-surface moisture – processes that interact intimately with both the macrophysical and microphysical components of the cloud scheme. Biases in low-cloud simulation or in the phase partitioning of mixed-phase clouds can therefore either amplify or dampen the simulated cloud loss and the downstream non-local Rossby wave responses.

Collectively, these parameterization uncertainties imply that the pronounced mid- and low-level cloud reduction and the resulting remote warming hotspots reported under idealized expansion scenarios should be interpreted cautiously as potentially lying at the upper end of the model-dependent range. To reduce these uncertainties, future work should prioritize targeted sensitivity experiments that systematically perturb key microphysical parameters (e.g., auto conversion thresholds, ice nucleation rates, and sub-grid variance formulations), employ higher horizontal resolution, or test alternative microphysics schemes (such as sectional or bin approaches). Improved observational constraints on cloud phase partitioning, liquid water path, and aerosol–cloud interactions over regions undergoing land-cover change will also be essential. In parallel, realistic assessments of future afforestation or bioenergy deployment require high-resolution, accurate land-cover datasets capable of reliably distinguishing forest from other vegetation types (Hoffmann et al., 2023). Generally, higher resolution means better performance in simulating cloud cover and key atmospheric dynamic processes (Kusunoki, 2026), leading to a more accurate understanding of the impacts of large-scale land-use change on cloud cover. In addition, dedicated checkerboard experiments or other validated separation frameworks could be employed in future studies to explicitly quantify the local and non-local components of LUCC-induced cloud responses (Chen et al., 2022; Zhao et al., 2026). Such experiments would complement the temperature-advection and wave-activity-flux diagnostics used here by providing a quantitative partitioning of cloud responses rather than only diagnosing the pathways of remote atmospheric adjustment. Taken together, our findings underscore that effective land-based climate mitigation strategies must extend beyond carbon sequestration targets alone. In addition to incorporating albedo-driven effects, as previous studies have emphasized (Wang et al., 2023), decision-making frameworks should also explicitly account for the potential climate risks arising from changes in cloud cover.

5 Conclusions

We use the current state-of-the-art fully coupled Community Earth System Model (CESM) to investigate the impact of two distinctly different future land use and cover changes (forest and bioenergy expansion) on long-term cloud cover trends. Simulation results show that future large-scale idealized forest and bioenergy expansions may exacerbate cloud-mediated regional climate warming, as they accelerate the loss of mid-low-level clouds while enhancing the increase of high-level clouds, thereby diminishing the overall climate-regulating benefits of cloud cover. Specifically, idealized afforestation leads to the fastest low-level cloud loss (1.14 times on global, 1.52 times on land), followed by bioenergy expansion (1.03 times on global and 1.23 times on land), whereas realistic afforestation mitigates such losses and generates positive climate effects by preserving or enhancing cloud cover. The loss of mid-low-level clouds amplifies the formation of regional warming hotspots through enhanced shortwave cloud radiative forcing trends, particularly in mid-to-high latitude regions of the Northern Hemisphere, where accelerated reductions in mid-low level cloud amounts align with amplified temperature trends (reaching up to 2.15 times the global average under idealized afforestation). The hemispheric asymmetry in cloud trends, quantified via Gini coefficients, further highlights how idealized forest and bioenergy expansions concentrate cloud loss trends within the latitudinal bands of those expansions, thereby potentially reshaping new climate warming hotspots. An interpretable machine-learning framework further identifies total precipitable water and relative humidity as influential predictors of the simulated cloud trends. Their high SHAP importance, together with the independent thermodynamic and boundary-layer diagnostics, is consistent with an important role of atmospheric moisture availability in the accelerated loss of mid- and low-level clouds under idealized forest and bioenergy expansion, which is also modulated by non-local factors. These results underscore that not all forest expansions yield climate cooling benefits; outcomes critically depend on the type of land conversion and the specific regions of implementation. This also highlights an important realization: even within bioenergy crops, the choice of crop function type is just as crucial as the choice of afforestation site. Therefore, policymakers should prioritize realistic afforestation strategies in suitable latitudinal zones to maximize carbon sequestration while avoiding unintended cloud-mediated warming.