Introduction

The incorporation of greater amounts of renewable energies into total national energy mixes is vital for limiting global warming to below 2 °C in the Paris Agreement1,2. Crucial in this transition is the extensive deployment of photovoltaic panels for low-carbon electricity production2,3,4,5. By the end of 2021, global solar PV capacity exceeded 800 GW, marking a tenfold increase over the past decade and contributing around 50% to the year’s renewable energy capacity growth6. Projections indicate that this capacity will exceed 18,200 GW by 2050 according to the International Renewable Energy Agency’s (IRENA) 1.5 °C Scenario7. This remarkable growth can be attributed to the reduced cost and improved efficiency of PV systems in harnessing carbon-free solar radiation and the growing commitment of nations to bolster their electricity portfolios with renewable energy sources. Despite PV’s important role in climate mitigation, a noteworthy concern is that the dark surface of solar panels, designed specifically to optimize solar radiation absorption, may lower terrestrial albedo at the PV sites and consequently result in a positive radiative forcing (RF), a warming effect that is at odds with efforts to curb anthropogenic warming through the reduction in greenhouse gas emissions8,9. Therefore, quantifying albedo changes and making clear the corresponding climatic consequences due to large-scale PV deployment is crucial for comprehending the complete role of PV in generating reliable forms of electricity that simultaneously address the issue of anthropogenic global warming.

However, understanding this potential drawback is hindered due to a lack of supporting radiation balance observations across a wide enough range of PV sites to quantify radiative feedback accurately. Previous studies have simulated the effects of PV deployment on climate10,11,12,13,14,15,16. Despite advancements in PV parmeterization10,11,12,13,14,15,16,17, many modeling studies12,13,14,15,16, when characterizing the PV’s effects on the surface energy budget, ideally assign overall terrestrial albedo values to regions featuring PV panel arrays, based on simplistic assumptions. While in-situ observation-based studies have provided valuable insights into changes in albedo due to PV construction18,19,20,21,22,23,24, such observations are confined to a limited number of specific local sites. Continuous satellite observations offer a unique opportunity to estimate terrestrial albedo changes by comparing albedo before and after solar farm construction or by contrasting solar farm albedo with that of adjacent regions25,26,27,28. Nevertheless, these investigations are still confined to limited PV sites. More importantly, the relatively smaller area covered by solar panels, compared with that of other landscape features, for example reservoirs29, poses challenges in accurately measuring the albedo change using mainstream optical products with spatial resolutions greater than the fine area covered by PV panels (Fig. 1b, c). Consequently, a comprehensive evaluation grounded on real-world observational data within a unified framework is imperative to assess the potential adverse albedo impact of photovoltaic farms.

Fig. 1: The spatial distribution of PV sites and conceptual illustration of linear regression method.
figure 1

a The spatial distribution of all PV sites (44,903) and filtered sites (352). b, c Depiction of a location featuring minimal PV panel coverage where MODIS grid albedo fails to characterize the albedo of the small PV-covered area. d Illustration of 4 PV polygons outlined by red solid lines, and the collection of grids with lower transparency represents the PV site. This background satellite imagery is sourced from ERSI World Imagery (2022-01-12) in ArcMap 10.8. e Demonstration of the correlation between land surface albedo and PV facilities’ area ratio across grids in the specific PV site depicted in d. The red solid line is the linear fitting line by the ordinary least-square method (OLS). The magnitude of the slope characterizes the albedo change caused by PV deployment (Eq. (5) in Methods). Besides, the upper and lower boundaries indicated by grey dashed lines, represent the site’s mean albedo values of background (with 0% coverage of PV facilities) and PV site (with 100% coverage of PV facilities), respectively.

Full size image

To overcome the abovementioned limitations and challenges, we introduce a linear parameterization method aimed at robustly detecting changes in land-surface albedo for individual PV farms, based on a global inventory of PV facilities30 and the Moderate Resolution Imaging Spectroradiometer (MODIS) albedo data product (Fig. 1). We first follow the methodology of Wohlfahrt et al.29 to calculate the blue-sky albedo of pixels at 500 × 500 m (MODIS grid cell; 134,520 cells in total) across 44,903 PV sites worldwide (Fig. 1a), and compare each site’s average albedo with its surroundings. However, direct comparisons entail substantial uncertainties in detecting PV-induced albedo changes because the coverage of PV facilities falls short of 100% for nearly all grid cells (134,297/134,520; Supplementary Fig. 1). Nonetheless, when sufficient grid cells exist within a PV site, a significant linear relationship between the PV-area ratio and albedo across these grid cells emerges (Fig. 1d, e); this relationship enables us to accurately quantify the albedo for both the background and the PV facilities at 352 representative sites globally (Fig. 1a). Next, we examine the factors influencing the albedo change, including Köppen–Geiger climate regime, background land cover and soil moisture. Finally, we calculate the albedo-induced global RF, which is further converted to anthropogenic carbon equivalence and compared with the negative RF equivalence from PV’s clean electricity generation (more details refer to Methods). Our findings reveal a PV-induced albedo decrease of −1.28 (−1.80, −0.90) × 10−2 (median and interquartile range), markedly smaller than assumed in global-scale climate simulations for widespread solar panel deployment in Earth system model (ESM) (refs. 14,15), supporting photovoltaic development for sustainable energy and climate goals.

Results

The overall decrease in terrestrial albedo caused by PV deployment

Across the 352 PV sites worldwide, the distributions of both background albedo and that following installation are significantly normal (Lilliefors test, P > 0.1; where P > 0.05 indicates significant normal distribution), with mean values of 0.1750 and 0.1606, respectively (Fig. 2a; Supplementary Fig. 2). The overall change in terrestrial albedo across PV sites is statistically significant but small (Fig. 2b; Paired t-test, P < 0.001), with a value of −1.28 (−1.80, −0.90) × 10−2, corresponding to a relative change of −7.0% (−10.8%, −9.0%) compared to the background (Fig. 2c, d). Of particular interest is the significant skewness observed in both absolute and relative albedo changes (Lilliefors test, P < 0.001). Skewing may be due to a bias in site selection for laying solar panels, with most sites (247/352) having a relative decrease of less than 10% in albedo.

Fig. 2: Analysis and comparison of mean albedo at PV sites and their corresponding backgrounds.
figure 2

a The frequency distribution of PV site’s mean albedo (with 100% coverage of PV facilities) and background albedo (with 0% coverage of PV facilities). µ and σ are the mean and standard deviation, respectively. b The boxplots of PV site albedo and background albedo. N denotes the number of sites. Paired t-test is used to test the significant difference between the PV site albedo with 100% coverage and background value (ns P > 0.05, *P < 0.05, **P < 0.01, ***P < 0.001). c, d show the distribution of albedo change magnitude (× 10−2) and percent change (%). Q25 and Q75 are 25th and 75th percent interval quantiles, respectively. The grey lines in c and d show the zero values. Lilliefors test is used to test the distribution’s normality, where a P-value larger than 0.05 suggests rejection of the null hypothesis that data comes from a distribution in the normal family at the 5% significance level.

Full size image

The albedo values of the PV sites, most of which exceed 0.1 (Fig. 2a), contrast with the simplistic assumed PV site albedo value (0.1) in simulations for global-scale climate feedbacks based on the ESM14,15 (Supplementary Table 1). Besides, the satellite-observed albedo changes are much smaller than those projected changes in ESMs14,15 (Fig. 2c, d; Supplementary Table 1). These disparities suggest that the assumptions in these modeling may inadequately represent the albedo and corresponding change at locations where PV panels are deployed. Furthermore, by comparing our findings with previous observation-based studies, we find general consistency in the direction of albedo change, while there are differences in magnitude (Supplementary Table 2). Differences between satellite observations and in-situ observations can be attributed to the settings of in-situ monitoring points for comparison19 and the spatial inhomogeneity of albedo within the grids of the MODIS albedo product31.

The changes in albedo following PV deployment vary across different regions (Supplementary Table 3; Supplementary Figs. 3 and 4). For example, the change in albedo at sites in the United States is −1.73 × 10−2, while sites in China exhibit a general albedo change of −1.23 × 10−2 (Supplementary Table 3). In India and other regions not in the United States or China, the reductions in albedo are relatively smaller, at magnitudes of −1.15 × 10−2 and −1.19 × 10−2, respectively. Albedo change exhibits relatively homogeneous patterns in longitude patterns, mainly ranging from −0.01 to −0.02, while the differences widen when latitude decreases from 45°N to 20°N (Supplementary Fig. 3; Pearson correlation, R = 0.20, P < 0.001). Notably, the most substantial albedo reduction (−9.19 × 10−2, −31.6%) happens in the United Arab Emirates (Supplementary Figs. 3 and 5). Moreover, there are a few sites with increased albedo (Fig. 2b, c; Supplementary Figs. 3 and 4), possibly due to high vegetation density or other anthropogenic causes of lower background albedo (Supplementary Fig. 6).

Impact of background environment on albedo change

Albedo changes significantly across different land-cover types after PV deployment (Fig. 3a–c; Paired t-test, P < 0.001), with distinct values, highlighting the land-cover type’s notable influence. The largest albedo decreases among all sites collectively is observed in sites located in open shrublands, having a value of −1.51 × 10−2 (Supplementary Table 3), followed by −1.42 × 10−2 for barren sites. Meanwhile, sites in croplands exhibit the smallest absolute decrease (−1.02 × 10−2), with a relative change value of −6.7%. Contrary to recent findings favoring greater albedo changes in barren areas (Supplementary Table 2; ref. 28), we uncover a larger overall albedo decrease in shrubland sites. This likely stems from barren sites being situated at higher latitudes (Supplementary Fig. 7), resulting in steeper solar panel angles, wider PV array spacing, and ultimately, a smaller fraction of the site covered by PV panels, leading to a reduced albedo change.

Fig. 3: The albedo change of sites covered by PV panels in different land-cover types and countries.
figure 3

a–c Boxplots of the background albedo (higher transparency) and the albedo in the site covered by PV panels (lower transparency) for different land-cover types with the paired points connected by gray line. The center line of each box represents the median; box limits represent the upper and lower quartiles; whiskers extend to 1.5 times the interquartile range. The captions show the median values of absolute albedo change (× 10−2). Gr, Cr, Ba, and OS represent sites in grasslands, croplands, barren and open shrublands, respectively. The numbers in parentheses after land-cover types represent the corresponding number of samples. Paired t-test is used to test the significant difference between the PV site’s mean albedo (with 100% coverage of PV facilities) and background albedo (with 0% coverage of PV facilities) (*P < 0.05, **P < 0.01, ***P < 0.001).

Full size image

Regionally, in China, which has the highest number of sites with land-cover types available for comparison, the maximum decrease in albedo is for barren areas (−1.27 × 10−2), followed by grasslands (−1.18 × 10−2) and croplands (−0.96 × 10−2) (Fig. 3b). In the United States, the PV-induced albedo changes exhibited the largest value in magnitude compared to other regions (−1.73 × 10−2), with larger decrease in albedo for grassland (−1.74 × 10−2), croplands (−2.55 × 10−2), barren (−2.58 × 10−2) and open shrublands (−1.70 × 10−2) compared to the same category in other regions (Supplementary Table 3). Moreover, we find that the categories of sites in different land-cover types from the United States notably differ from those of China and India (Supplementary Table 3). This indicates that, despite consistent land-cover types, the reduction in albedo at PV sites exhibits notable spatial variation, suggesting the influence of factors beyond land-cover types.

Further analysis reveals that climate regime plays a pivotal role in influencing albedo changes, even when considering the same land cover type (Fig. 4 and Supplementary Fig. 8). For example, sites located in barren have significantly different albedo reductions between the BWh (desert, arid and hot; −2.57 × 10−2) and the BWk climate classes (desert, arid and cold; −1.24 × 10−2) (Wilcoxon rank sum test, P < 0.001). Additionally, sites located on water bodies in cold winter climates (Dwb) exhibit a higher absolute albedo decrease (−2.92 × 10−2; −23.9% of the background) compared with those in the Cfb regime (temperate, no dry season, warm summer). This effect may be due to solar panels reducing the high albedo of freezing water during the cold season. We also explore whether different climates cause varying albedo changes across countries. A comparison of sites over grasslands in China and the United States, where sufficient samples are available, reveals that nearly 25% of PV sites over grasslands in China are located under cold and dry winter conditions (Dwa, Dwb, and Dwc regimes; Supplementary Fig. 9a, b), with a median albedo change of −1.02 × 10−2. In contrast, no such sites exist in the United States, potentially contributing to the lower albedo change at PV sites over grasslands in China. Nonetheless, even for grassland sites under similar climatic conditions (e.g., BSk regime), the albedo changes at photovoltaic sites in the two countries differ (−1.27 × 10−2 in China; −1.74 × 10−2 in the United States). This disparity could be attributed to the different PV array spacing induced by variations in latitude (Supplementary Fig. 9c).

Fig. 4: Analysis of absolute albedo change and relative albedo change at PV sites across different land cover types and climate regimes.
figure 4

a The ratio of site numbers in specific climate regimes to the total sites of corresponding land cover types. WB, Ba, Cr, OS, Gr, Sa, and Wsa represent sites with the background land cover of water bodies, barren, croplands, open shrublands, grasslands, savannas, and woody savannas, respectively. Additionally, the nomenclature following each land cover type signifies the climate regimes, which are derived from the Köppen–Geiger classification map, and more details are shown in Supplementary Table 5. Only groups with a site number exceeding three are depicted in this figure. b, c The box plots of albedo change (×10−2) and relative albedo change compared to the background albedo, respectively. The white center line of each box represents the median; box limits represent the upper and lower quartiles; whiskers extend to 1.5 times the interquartile range. Captions near the box plot show the median and interquartile range (IQR) of corresponding categories, respectively.

Full size image

Other factors could also influence PV-induced albedo change, including soil moisture. We find that site-averaged soil water content in our samples averaged approximately 0.1 (Supplementary Fig. 10), showing a significant negative (Pearson) correlation with both initial and altered albedo (P < 0.001; Supplementary Fig. 11a). This wetness-albedo correlation aligns with prior studies32,33, where increased soil moisture darkens surfaces due to geometric effects and enhanced forward scattering34,35. Meanwhile, higher soil moisture would lead to greater vegetation biomass that absorbs more incoming radiation, which further lowers albedo. However, the linear response of albedo change to soil moisture content is not statistically significant (P > 0.05; Supplementary Fig. 11b). This observation holds true in most cases, even when considering sites in specific climate zones with corresponding land cover types (Supplementary Fig. 12), indicating that the potential impact of soil water content on albedo change may be minimal.

The radiative forcing and carbon equivalence due to albedo change

The reductions in terrestrial surface albedo following large-scale solar panel deployment will induce a positive global RF at the top of the atmosphere, consequently, a warming effect that counteracts the negative radiative forcing due to the substitution of clean solar generation for fossil fuels8. Overall, the total global RF due to albedo change across 352 sites is approximately 7.48 μW m−2, which is equivalent to anthropogenic carbon emissions of 2.77 Tg C (Fig. 5a). Individually, the values range from −0.03 µW m−2 to 0.52 µW m−2 (Supplementary Fig. 13). The majority (345/352) exhibits positive values, compared with only 7 sites showing negative values indicating potential cooling. Among those sites with positive RF, we observe large variability in the values, spanning several orders of magnitude. For instance, the RF values at sites in the United States are notably higher than those at sites in Europe (Fig. 5a).

Fig. 5: The global radiative forcing (RF) and carbon equivalence (CE) due to albedo change.
figure 5

a The spatial pattern of the global RF caused by PV deployment. The insert shows the top 30 sites’ RF values alongside corresponding anthropogenic carbon equivalence. b–d The relationship between three key variables—albedo change, mean downward shortwave radiation, PV site area—and the global RF at the top of the atmosphere. The relative differences are expressed as the absolute percentage changes of each variable relative to its respective minimum absolute value. The black scatters show the relationship between RF and relative difference of corresponding variable, while the upper bars represent the frequency distribution of relative difference. The captions show the Pearson partial correlation coefficients between RF and each variable (*P < 0.05, **P < 0.01, ***P < 0.001), respectively. 7 sites with negative RF are not shown in the figure above.

Full size image

Global RF is determined by three main factors: the area covered by PV panels, the albedo change caused by PV deployment, and the mean downward shortwave radiation in the region. The area’s impact on RF is more pronounced due to its extensive variability across multiple orders of magnitude, compared to albedo change and radiation (Fig. 5b, c). Furthermore, there is an interaction between the three variables (Supplementary Fig. 14), indicating that larger PV farms in sunnier locations, are often accompanied by larger albedo reductions, thereby magnifying the resulting global positive RF. We further examine the local RF (Supplementary Methods), crucial for regional energy budget, which ranges from −4.48 W m−2 to 20.56 W m−2 (Supplementary Fig. 15). Notably, the desert site in the United Arab Emirates exhibits the largest positive local RF value, because of its exceptionally large albedo change compared to other sites (Supplementary Figs. 5, 15 and 16), suggesting that deploying PV on desert land could lead to a larger temperature disturbance.

Discussion

Understanding the unintended climate impacts of widespread solar panel deployment is crucial for tackling climate change. It’s essential to practically characterize the albedo changes at PV sites to refine climate models at both global and regional levels. The overall surface mixed albedo of a PV farm reflects both the reflectivity of solar panels and that of the natural surface, accounting for the required spacing between arrays (ref. 36; Supplementary Fig. 17). Given that the albedo of most land cover types exceeds 0.1 (ref. 37), neglecting the background albedo in spacing in some global-scale ESM-based simulations14,15, which utilized simplified fixed albedo of 0.1 to represent PV sites over the desert, can lead to lower mixed albedo at PV sites compared to observations (Supplementary Tables 1, 2 and 4). This, in turn, results in a larger relative albedo change from the background (up to 75% decrease; refs. 14,15), and thus an overestimated climate response.

Introducing the packing factor, a parameter representing the percentage of interested land covered by panels in the PV site10,38, into global-scale Earth system simulations offers a straightforward method to address this concern. It enables the mixed albedo of regions with PV installations to more accurately align with observed values and refine the heterogeneity in PV-induced albedo change caused by the underlying background characteristics. Moreover, while regional climate modeling allows for more complicated PV parameterizations10,11,12,13, some schemes characterize the site’s mixed albedo based on only a few single-point field measurements, introducing scale-related uncertainties (e.g., ref. 12). Our analysis provides valuable insights for future studies on climate feedbacks associated with PV installations across diverse land covers under different climate conditions, despite focusing on the mean state change of albedo caused by PV deployment and ignoring the variable albedo of PV site (Supplementary Table 4).

The overwhelmingly widespread deployment of solar panels in the future will further change the Earth system’s energy budget and could result in larger positive global RF. The global effect of radiative forcing caused by the 352 selected PV sites, which represent approximately 24% of the total area (44,903 sites), is approximately 7.48 µW m−2 (Methods). Assuming that magnitudes of other unselected PV sites’ radiation and albedo change remain consistent with the selected ones, the total global RF of all sites (44,903) could be around 30.79 µW m−2. By 2050, according to the projected installed solar PV capacity of exceeding 18,200 GW (~37 fold the capacity in 2018) in the IRENA’s 1.5 °C Scenario7, the global RF would potentially reach more than 1135 µW m−2 (equivalent to anthropogenic carbon emissions of approximately 426 Tg C), compared with 3300 µW m−2 obtained from the idealized assessment under a similar scenario of PV installation capacity8.

The positive global RF due to the reduction in albedo impairs the emission reductions of solar panels and increases the energy payback time in the life cycle9. Therefore, understanding how this adverse albedo effect counteracts the emission reductions from solar generation is of paramount importance. Annual generation at PV sites varies from 2.84 × 107 to 4.74 × 109 kWh year−1, while electrical output per unit area ranges from 70.06 to 79.94 kWh year−1 m−2 (Supplementary Fig. 18). Despite these variations, the clean electricity generated by PV in most sites could offset their adverse albedo impacts within a single year (break-even time; Supplementary Fig. 19). This indicates a cooling effect in the subsequent years of PV operation, emphasizing the positive role of deploying PV panels in mitigating global warming. Nonetheless, our estimation of break-even time is idealized and does not include several specific factors that could potentially prolong this period. These factors include the omission of other PV-related radiative forcing, such as longwave forcing, and the use of idealized PV generation calculations involving overlooking the degradation of PV generation efficiency over time. Additionally, we do not consider the carbon sequestration changes in natural lands caused by PV installations, as these are relatively minor compared to the carbon offsets at the PV site (Supplementary Fig. 20). Furthermore, the short-term energy transformation from solar radiation to electricity on local temperature disturbance and broader global temperature responses are not comprehensively explored here. Clarifying these complexities could be a potential avenue for future research.

Strategic planning is crucial to the climatic benefits of future large-scale PV deployment. Transitioning lands to PV farms requires optimizing PV generation per unit area and minimizing the albedo reduction to shorten break-even times. Utilizing more efficient solar panels increases electrical output per area and land-use efficiency39, thereby reducing the break-even time through enhanced carbon avoidance (Supplementary Methods) and decreased positive global RF due to smaller land requirements. Moreover, integrating PV on low-albedo land-cover types, such as croplands and grasslands, not only mitigates albedo effects and optimizes land use but also offers additional benefits. For instance, in some dryland agrivoltaic systems, crops benefit from reduced land-surface temperatures and optimized soil water balance, leading to higher yields, while solar panels sustain elevated conversion efficiencies40,41. Additionally, installing solar panels on water surfaces in regions without freezing concerns brings year-round electricity generation benefits. Floating PV systems combined with traditional and pumped storage hydropower not only generate substantial electricity but also curb water evaporation42.

However, the deployment of PV panels also carries potential environmental and ecological risks24. Changes in carbon sequestration from PV installations on natural lands, though might be minor compared to the carbon avoidance of generation, are unneglectable compared to the land’s original state (Supplementary Fig. 20a, c). This is mainly due to landscape reshaping38, influencing local native vegetation dynamics and soil microbial characteristics43. Consequently, ecologically rich lands and vital ecosystems should be avoided by the energy industry44. Additionally, in certain croplands requiring high solar radiation or day-night temperature difference, the shading of solar panels reduces crop yield and quality45,46. Floating PV systems may also influence water quality42, warranting comprehensive impact studies. In relative terms, converting highly degraded barren to a solar farm, despite suffering from its positive radiative forcing and potential extension of energy payback time, may be more cost-effective when considering land and ecosystem service values, making it a suitable priority target for conversion. Therefore, future PV expansion requires careful consideration to maximize the climatic benefits and minimize ecological disruptions and environmental influences.

Methods

The datasets of PV facilities, albedo, and underlying surface information

Our study utilizes a dataset of global PV facilities derived from global Sentinel-2 and SPOT6/7 imageries (from 2016-06-01 to 2018-09-30; ref. 30). This PV dataset comprises 68,661 vector polygons delineating the solar generating units’ coverage worldwide, along with corresponding information containing estimated capacity and area. After eliminating duplicates, 68,655 polygons remained. Furthermore, we supplemented the missing attributes of country information for each polygon by matching them with neighboring polygons containing available country data.

To align with the temporal coverage of imageries used in generating the PV dataset and improve the number of available information, we downloaded the Breathing Earth System Simulator (BESS) terrestrial shortwave radiation47 and the MODIS MCD43A1 Version 6.1 BRDF/Albedo Model Parameter dataset48 spanning from 2019 to 2021. These two datasets, offering daily information at spatial resolutions of 5 km and 500 m, respectively, were used to calculate the blue-sky albedo of grid cells at 500 m resolution (Eqs. (1)–(4)).

Additionally, we conducted a comprehensive assessment of the land cover and climate regimes in the locations where the extracted samples are situated. As for land-cover types, we used the International Geosphere-Biosphere Programme (IGBP) classification in 2019 derived from the MODIS Land Cover Type Product (MCD12Q1 V006) (ref. 49). This product, available at a spatial resolution of 500 m and with an annual time step, offers detailed classifications of land cover (17 types). Besides, we identified the climate zones based on the global map of the Köppen–Geiger climate classification (1980–2016) at a 1-km resolution50, with more detailed information available in Supplementary Table 5.

We further used soil moisture data to explore the effect of environmental dryness and wetness on albedo change. For soil moisture information, a global map of downscaled Soil Moisture Active Passive (SMAP) dataset at a 1-km spatial resolution was used in our analysis51. This product offers daily soil moisture information at two specific time points, namely 6:00 a.m. and 6:00 p.m. It demonstrates high reliability in mid-and-low latitude regions, aligning well with the spatial distribution of PV facilities. To match the time scale of other data used in our study, daily soil moisture data was downloaded spanning from 2019 to 2021 and averaged over the 3 years.

The calculation of actual albedo

The calculation of actual albedo, or blue-sky albedo at a specific grid (500 × 500 m) involves a straightforward combination of white-sky albedo and black-sky albedo, which can be written as follows52:

$${{Albedo}}_{{mix}}={{Albedo}}_{{mix}{\_}{ws}}\gamma +{{Albedo}}_{{mix}{\_}{bs}}\left(1-\gamma \right)$$
(1)

where \({{Albedo}}_{{mix}}\) is the satellite-based terrestrial albedo at grids with mixed surface features at hourly time steps, \({{Albedo}}_{{mix\_ws}}\) and \({{Albedo}}_{{mix\_bs}}\) signify the white-sky and black-sky albedo at hourly time steps, respectively. The variable \(\gamma\) is the ratio of the terrestrial diffuse shortwave radiation to the terrestrial total downward shortwave radiation. This ratio is an empirical function of the solar zenith angle53, which can be expressed as:

$$\gamma =\left\{\begin{array}{cc}0.943+0.734\rho -4.9{\rho }^{2}+1.796{\rho }^{3}+2.058{\rho }^{4} & \rho \, < \, 0.8\\ 0.13 & \rho \, > \, 0.8\end{array}\right.$$
(2)

where \(\rho\) is a ratio of hourly global terrestrial shortwave radiation (W m−2) to the product of the solar constant (1367 W m−2) and \(\cos {\theta }_{i}\), where \({\theta }_{i}\) is the solar zenith angle at a certain hour (\(i\)). The hourly land-surface shortwave radiation values were derived by scaling the daily average radiation, as provided by the BESS radiation product47, against the daily average extraterrestrial radiation. This ratio adjusts the daily radiation values to an hourly scale, reflecting variations in extraterrestrial radiation throughout the day, under the assumption of consistent atmospheric conditions29.

The Bidirectional Reflectance Distribution Function (BRDF) characterizes the reflective properties of an object’s surface, determining how incident light is scattered and distributed in different directions upon reflection. The satellite-based albedo product MODIS MCD43A1 V061 provides three daily BRDF albedo parameters for shortwave broadbands48, which can be used to calculate the white-sky albedo and black-sky albedo of the target pixel by using the following equations, respectively:

$${{Albedo}}_{{mix}{\_}{ws}}=f_{{iso}}+0.189184f_{{vol}}-1.377622f_{{geo}}$$
(3)
$$\begin{array}{c}{{Albedo}}_{{mix}{\_}{bs}}=f_{{iso}}\\ +\, \left(-0.007574-0.070987{\theta }_{i}^{2}+0.307588{\theta }_{i}^{3}\right)f_{{vol}}\\ +\, \left(-1.284909-0.166314{\theta }_{i}^{2}+0.041840{\theta }_{i}^{3}\right)f_{{geo}}\end{array}$$
(4)

where for both equations, \(f_{{iso}}\), \(f_{{vol}}\) and \(f_{{geo}}\) are three corresponding daily BRDF parameters of the target pixel derived from the albedo product.

Here, we used a QC value of 0 according to the QC-flag layer provided by the MODIS albedo product, as well as a set threshold of 70 degrees for solar zenith angle54, to get pixels with high-quality BRDF parameters. We then followed Wohlfahrt et al.29 to calculate the hourly blue-sky albedos for grids containing PV facilities. The hourly grid albedo values were subsequently used to derive daily, monthly and 3-year (2019–2021) weighted averages by utilizing corresponding time-scale downward shortwave radiation.

Albedo change due to PV panel deployment

We created 44,903 PV domains based on 68,661 PV vector polygons, symbolizing different sites (Fig. 1a and Supplementary Fig. 21). The PV vector polygons were rasterized based on the grid cell of 500 m, mirroring the spatial resolution of MODIS albedo data (Fig. 1d and Supplementary Fig. 21). The resulting isolated PV domain is a collection of grids containing PV facilities that is representative of the region with the highest number of PV polygons. In our preliminary analysis, we considered the average albedo within a buffer area of 2 grids around each domain as the background albedo (Supplementary Fig. 21). Consequently, the mean albedo within the PV domain represented the altered albedo caused by the replacement of vegetation with PV panels. Problematic is that this method tends to introduce large uncertainties because the size of individual PV polygons typically does not cover an entire grid.

Therefore, we hypothesized a linear relationship between the area coverage of PV facilities range within a grid and the corresponding grid albedo at a site. We proceeded under the assumption that the random errors within this relationship can be ignored, allowing us to conduct an ordinary least-squares linear fit (Fig. 1e). The result indicates a significant relationship (P < 0.05) when the number of pixels in a PV site is large enough. The slope and intercept from the regression parameters indicate albedo change and background albedo, respectively. Hence the following equation applies:

$${\overline{{Albedo}}}_{{mix}}=\left({\overline{{Albedo}}}_{{PV}}-{\overline{{Albedo}}}_{{background}}\right){AR}+{\overline{{Albedo}}}_{{background}}$$
(5)

where \({AR}\) represents the PV facilities’ area ratio of the grid, which was calculated by using the ArcPy package in Python 2.7. \({\overline{{Albedo}}}_{{mix}}\) is the 3-year (2019–2021) weighted albedo average of a grid with mixed features, \({\overline{{Albedo}}}_{{PV}}\) and \({\overline{{Albedo}}}_{{background}}\) represent the albedo at grids completely filled by PV facilities (\({AR}\) = 1) and by background land cover surface (\({AR}\) = 0), respectively. Furthermore, we compared albedo change values calculated by direct comparison (PV site albedo – buffer albedo) and by this linear parameterization method (Supplementary Fig. 22). Our analysis showed that ~82% of sites (289/352) exhibited a positive relative difference (with a median value of 37%) between albedo change values calculated with and without PV fractions, suggesting a notable underestimation in the albedo change calculated by direct comparison.

The collection of representative sites covered by solar panels

A two-step filtering process was performed to ensure the reliability of independent parameters for each site, due to either (a) insufficient pixels for analysis in many sites or (b) clustering of area ratios of different PV pixels around a certain value in some sites, leading to unreliable regression parameters. The initial filtering criteria comprised the following: (1) the number of pixels in a domain should exceed 10; (2) the difference between the maximum and minimum area ratio values across all pixels within an individual PV site should be larger than 0.5; (3) the P value of F-test for the linear model is required to be smaller than 0.05. To evaluate the sensitivity of the number of filtered sites to the criteria, we modified the criteria to observe the resulting changes. The number of filtered sites is sensitive to the setting of increasing the grid cell’s number in a site, but shows less evident changes concerning the difference between maximum and minimum area ratios of PV facilities’ range to the corresponding pixel is not obvious (Supplementary Table 6). Therefore, we finally chose the above criteria to acquire sites with robust slope and intercept parameters. After the first filtering, 412 sites remained eligible for subsequent analysis.

Further, we re-examined the filtered samples to ensure the precise detection of the PV vector polygons’ detection and their corresponding boundaries. We utilized the satellite imagery dated 12 January 2022 from ESRI World Imagery in ArcMap 10.8 to conduct visual interpretation. In cases where regions in this image were unclear, additional detection was carried out using Google Earth. Samples with incorrect boundaries but minimal impact on the regression analysis were included. In the end, 352 representative sites (approximately 1174 km2) accounting for ~24% of the total area covered by PV facilities were retained for further analysis. The limitations posed by the coarse-resolution radiometric products restrict a more extensive investigation into albedo changes at additional sites. However, researchers have already tried to explore the inversion of higher spatial resolution albedo products55,56,57. The general applicability of our approach allows for easy adaptation when such products become available publicly.

Overall, most PV sites, both in the original PV sites (99.0%) and in post-selected sites (92.9%), are located in mid-low latitude regions of the northern hemisphere (from 0° to 60°N). We conduct aggregate statistics on our samples at both the country and land cover levels. While our samples span 26 countries, the distribution is uneven. Notably, China hosts the highest number of sites, accounting for over half of the total (188 sites; Supplementary Table 3), with the majority (95.7%) located in the northern part of the country. A substantial number of sites are also found in the United States (66 sites, 18.8%) and India (30 sites, 8.5%). Regarding land cover, our samples encompass regions demarcated as grasslands (144 sites), croplands (61 sites), barren areas (77 sites), open shrublands (39 sites), and other types (31 sites) (Supplementary Tables 3 and 7).

Global radiative forcing and carbon equivalence due to albedo change

Global radiative forcing (\({RF}\); W m−2) at the top of the atmosphere reveals the change in the Earth system’s energy balance8. It is defined as the difference between the sun’s incoming energy into the Earth system and the energy radiated back into space. Positive RF may imply a warming effect on the entire Earth’s surface. A typical simplified model is applied here to calculate radiative forcing58,59. We assume that the global effect of PV RF due to albedo change is instantaneous29. Nevertheless, the characterization of instantaneous RF relies on the mean albedo change (2019–2021), derived through the linear parameterization method based on the 3-year weighted grid albedo values (2019–2021). Hence, the RF of radiance imbalance from albedo change can be quantified as follows:

$${{RF}}_{\overline{\Delta {Albedo}}}=\frac{{\bar{R}}_{{SR}}^{\downarrow }\overline{\Delta {Albedo}}{T}_{{SR}}^{\uparrow }{A}_{{PV}}}{{A}_{E}}$$
(6)

where \({\bar{R}}_{{SR}}^{\downarrow }\) is the 3-year average incident shortwave radiation (2019–2021) at the terrestrial surface (W m−2), \(\overline{\Delta {Albedo}}\) is the mean albedo change due to PV deployment, which is calculated from the 3-year weighted average grid albedo (2019–2021) by using the linear parameterization method (Fig. 1e), \({A}_{{PV}}\) represents the scope area covered by PV facilities in a PV site, \({A}_{E}\) denotes the Earth’s surface area (510 × 106 km2), and \({T}_{{SR}}^{\uparrow }\) is the upward transmittance constant, set at 0.854 (ref. 60).

To quantify the anthropogenic carbon emission equivalence of global RF, we introduced the equation of Myhre et al.61:

$${{RF}}_{\overline{\Delta {Albedo}}}={RE}\cdot {{ln}}\left(1+\frac{\Delta C}{{C}_{0}}\right)$$
(7)

where RE is the radiative efficiency (5.35 W m−2), \(\Delta C\) is the change of CO2 concentration (ppm) in the atmosphere based on reference concentration (410 ppm for this study). Assuming a small \(\frac{\Delta C}{{C}_{0}}\), Taylor series expansion is applied to first-order linearization of Eq. (7):

$${{RF}}_{\overline{\Delta {Albedo}}}={RE}\frac{\Delta C}{{C}_{0}}$$
(8)

And introducing the airborne fraction (\(\rho\) = 0.44, ref. 60), then the carbon equivalence (\({CE}\); g C) of albedo-induced radiative forcing can be calculated as follows:

$${{CE}}_{\overline{\Delta {Albedo}}}=\frac{\Delta C}{\rho }=\frac{{{kC}}_{0}{{RF}}_{\overline{\Delta {Albedo}}}}{\rho {RE}}$$
(9)

where \(k\) is the constant factor of converting ppm to g C (2.13 × 1015).

However, part of the incoming radiation is converted into electricity by solar panels. Here we assumed this was a carbon-free process and ignored any degradation in generation efficiency. Therefore, the reduced carbon emissions each year from each PV site’s solar generation compared to coal-fired plants can be derived as follows:

$${{CE}}_{{{\rm{gen}}}}=8760\lambda \cdot {CI}\cdot {CF}\cdot {Cap}$$
(10)

where 8760 is the corresponding number of hours in a year, \(\lambda\) is the coefficient (0.27) transforming the carbon dioxide intensity (g CO2 kWh−1) to carbon intensity (g C kWh−1). \({CI}\) is the carbon dioxide intensity (900 g CO2 kWh−1) of coal-fired plants in 2018 (ref. 62), \({CF}\) is the mean capacity factor (0.11) of solar PV in the world63 and \({Cap}\) (kW) is the total capacity of a PV site, which is the sum of estimated nominal peak alternating current generating capacities of each solar generating units in the site. Each solar generating unit corresponds to a vector polygon in the global PV dataset, where the capacity of each unit has been evaluated based on its size, the efficiency of the solar panels, and other factors30.

Then the number of years needed to offset the global RF caused by albedo change is derived by combining Eqs. (9) and (10):

$${\mbox{`}} {Break}\; {even}\; {time}{\mbox{'}} =\frac{{{CE}}_{\overline{\Delta {Albedo}}}}{{{CE}}_{{{\rm{gen}}}}}$$
(11)