Albedo changes caused by future urbanization contribute to global warming

Land cover

Urban landscapes are characterized by small clusters of patches, forming land mosaics that are distinct from natural landscapes. An accurate estimation of climate forcing requires a land cover dataset at high resolutions that does not omit small urban patches. In this study, the RF estimates are based on 500-m and 1-km land cover datasets. This fine resolution is necessary to preserve spatial details of small urban patches while avoiding the large underestimation of urban land areas at coarse resolution (e.g., ~19% underestimation at 10 km compared to that at 1 km)³. We used 500-m resolution MODIS Land Cover product (MCD12Q1v006) for historical land cover changes. For future urban land cover distributions, we used the global urban land expansion products simulated under the SSPs for 2030–2100 (i.e., Chen-2020)⁴. The simulation performance was tested using historical urban expansion from 2000 to 2015 based on Global Human Settlement Layer⁵¹, where the agreement between simulated and observed urban land was evaluated using the Figure of Merit (FoM) indicator⁵² that has showed similar or better values than those reported in other existing land simulation applications⁴. The high-resolution Chen-2020 also shows very high spatial consistency with the prominent coarse resolution global urban land projection LUH2 that is recommended in CMIP6⁴. Considering different scenarios is also necessary to account for the uncertainties of future socioeconomic and environmental conditions, so we included simulated urban lands under three scenarios (Supplementary Table 1): Sustainability -SSP1, Middle of the Road – SSP2, and Fossil-fueled Development – SSP5⁵³. Within each SSP scenario, the product provides a likelihood map of each grid becoming urban, based on 100 urbanization simulations. We used the likelihood map to account for spatial uncertainties of urban expansion by deriving 90% confidence intervals of projected urban land demand within a SSP scenario. We used the MODIS IGBP Land Cover classes (Supplementary Table 2) and resampled the original 500-m resolution MODIS products in 2018 to 1-km resolution to match the future simulations when it was used as a baseline year. To isolate the independent effect of urbanization (vs other types of land uses) in future estimates, land covers that are not converted to urban are assumed to have the same cover types as in 2018 (i.e., the baseline year). Though there are other global land cover products for current periods, we choose the MODIS IGBP land cover products because the albedo look-up maps (LUMs) were based on IGBP land cover types (see Albedo Look-Up Maps).

To further evaluate the uncertainties caused by different projections of future urbanization, we also included the other two SSPs from Chen-2020, and another two 1-km resolution urban land cover products projected for the future for the purpose of comparison. The other two products include four projections of SRES scenarios (i.e., A1, B1, A1B, and B2) (i.e., Li-2017 mentioned above)³ and one without scenario description but assumed historical development would continue (i.e., Zhou-2019 mentioned above)². These projections of future urban land expansion were calibrated with different historical urban land products and can be regarded as independent.

Albedo look-up maps (LUMs)

Albedo Look-Up Maps (LUMs)³¹ were derived from the intersection of MODIS land cover⁵⁴ and surface albedo⁵⁵ products, which are used to determine the albedo values for each IGBP land cover type by month and by location. Monthly means of white-sky (i.e., diffuse surface illumination condition) and black sky (i.e., direct surface illumination condition) during 2001–2011 were processed for snow-free and snow-covered periods for each of the 17 IGBP land cover classes at spatial resolutions of 0.05°−1°³¹. The LUMs have been verified by comparing the reconstructed albedo using the LUMs with the original MODIS albedo, which shows very similar values³¹. We used the LUMs at a resolution of 1° due to the significantly fewer missing values, to assure the spatial continuity of albedo changes at a global scale while keeping the matches with the 1° resolution of radiation data and RF kernels. The underlying assumption is that albedo of the same land cover type varies insignificantly within a 1° grid.

Snow and radiation product

Snow cover can significantly change the albedo of land regardless of cover types (Supplementary Fig. 4). In this study, we tally monthly albedo using snow-free and snow-covered categories in estimating RF. Past and present snow-free and snow-covered conditions were derived from level 3 MODIS/Terra Snow Cover (MOD10CM.006)⁵⁶ at 0.05° spatial resolution and resampled to a 1° spatial resolution. Monthly means of 2001–2005 vs 2015–2019 were used for 2001 and 2018 respectively. For future periods, ensemble mean snow cover for each year and month, projected under the CMIP5 framework for three Representative Concentration Pathway (RCP) scenarios (i.e., RCP2.6, RCP4.5, and RCP8.5) were used (for more details see Supplementary Note 2B). By comparing the model outputs with MODIS observations for a recent decade (2006–2015), we found that the multi-model mean snow cover was systematically biased compared to MODIS observations. Consequently, we calibrated the ensemble mean projections by subtracting the biases for the grids. In each 10^th year of the future (e.g., 2030, 2040, etc.), the decadal monthly mean snow cover (e.g., 2026–2035 for 2030, and 2036–2045 for 2040, etc.) was used for the year.

We used the long-term monthly averages (1981–2010) of diffuse and direct incoming surface solar radiation reanalysis Gaussian grid product from National Centers for Environmental Prediction (NCEP)⁵⁷. Visible and near infrared beam downward radiation and diffuse downward radiation from NCEP were used to compute the white-sky and black-sky fractions. As for snow cover, ensemble mean shortwave radiation at surface (SW^SF) and at top-of-atmosphere (SW^TOA) projected from CMIP5 models (Supplementary Note 3C) for RCP2.6, RCP4.5, and RCP8.5 were collected for empirically computing future albedo kernels (see section 3.4 below).

Radiative kernels

Radiative kernels were used to compute top-of-atmosphere RF due to small perturbations of temperature, water vapor, and albedo. We used the latest state-of-the-art albedo kernels calculated with CESM v1.1.2⁵⁸ to compute RF in 2018 relative to 2001. In brief, the albedo kernel is the change in top-of-atmosphere radiative flux for a 0.01 change in surface albedo. The CESM1.1.2 kernels are separated into clear- and all-sky illumination conditions. We used the all-sky kernels because we include both black-sky and white-sky albedos. For future periods, because there are no available radiative kernels produced from general circulation models, we approximated the future kernels using an empirical parameterization following Bright et al.⁵⁹:

where m is the month, i is the location, and SW^SF and SW^TOA are the surface and top-of-atmosphere shortwave radiation; dividing by −100 is for matching the CESM1.1.2 kernel definition of a 0.01 change in surface albedo.

Estimation of albedo change and RF

We analyzed the RF in 2018 due to albedo changes caused by urbanization since 2001 (2018–2001), and in the future from 2030 to 2100 at decadal intervals (i.e., 2030, 2040, 2050, …, and 2100) since 2018 under three illustrative scenarios: SSP1-2.6, SSP2-4.5, and SSP5-8.5, which combine SSP-based urbanization projections and RCP-based climate projections. The three illustrative scenarios were selected following the scenario designation of the latest IPCC report⁵⁰ and represent low greenhouse gas (GHG) emissions with CO₂ emissions declining to net zero around or after 2050, intermediate GHG emissions with CO₂ emissions remaining around current levels until the mid-century, and very high CO₂ emissions that roughly double from current levels by 2050, respectively. We selected 2018 as the baseline year to divide the past from the future because 2018 was the latest year with available MODIS land cover products at the time of this study. We used ArcGIS 10.6 to produce spatial maps of all variables, including area of each land cover type within a 1° × 1°-grid, snow cover, albedo, radiation, and kernels, and R 3.6.1 to compute the RF.

We focused only on albedo changes induced by urbanization, including the conversions from all other 16 IGBP land cover types to urban land. The changes of albedo for each grid (x, y) of a month (m) were obtained by computing the difference between albedo of that grid in the baseline year (t = t₀) and in a later year (t = t₁) with urban expansion:

where α_{m, t = t1} (x, y) and α_{m, t = t0}) (x, y) is the albedo for each grid (x,y) of a month (m) at the base year and later year respectively; the grid-scale albedo is computed as the weighted sum of albedo by land cover types with the weighing factor corresponding to areal percentage of a land cover within the grid. The albedo for each land cover type of a grid was then obtained by applying the albedo LUMs that provide spatially continuous black-sky, white-sky, snow-covered, and snow-free albedo maps for a given month for each land cover. Firstly, monthly mean albedo is computed as:

where m is the month, t is the year, l is the land cover type, f_l is the proportion of a cover type within the grid, f_s,m,t is the fraction for snow-covered (s = 0) and snow-free (s = 1) conditions of the time (m, t), f_r,m,t (x, y) is the fraction for white-sky (r = 0) or black-sky (r = 1) conditions of the time, and α_l,s,r,m (x, y) is the albedo for land cover type l in month m that is extracted from the albedo LUMs corresponding to snow condition (s) and radiation condition (r). The annual mean albedo change is reported as the mean of monthly albedo change:

The conversion of other land covers to urban land can contribute differently to the global RF, as the total area that is converted into urban land is different among non-urban land covers and the albedo differences between urban land and non-urban land cover types vary. To estimate the proportional contributions of different land conversions, we first decomposed the total albedo of each grid into the proportion of each land cover type:

The global RF due to albedo change caused by conversion from each non-urban land cover type (l ≠ 13) to urban land (l = 13) (see Supplementary Table 2 land cover labels) was calculated as:

where i refers to a grid, n is the total number of pixels on global lands, A_Earth is the global surface area (5.1  ×  10⁸ km²), α_{13,m,t = t1}) (i) is the albedo of urban land in month m in the later year with urban expansion, α_{l,m,t = t0} (i) is the albedo of a targeted non-urban land cover type in the base year t₀, Δp_l→13 is the percentage of the non-urban land cover type that is converted to urban land in the year t₁ compared to year t₀, Area(i) is the area of the pixel, and K_m (i) is the radiative kernel at the grid.

The global RF due to urbanization-induced albedo changes was then calculated as:

GWP: CO₂-equivalent

We followed GWP calculations by explicitly accounting for the lifetime and dynamic behavior of CO₂ to convert RF to CO₂ equivalent^60,61:

where k_CO2 is radiative efficiency of CO₂ in the atmosphere (W/m²/kg) at a constant background concentration of 389 ppmv, which is taken as 1.76  ×  10¹⁵ W/m²/kg⁶², and RF_∆α,global is the global RF caused by albedo changes (W/m²). ({y}_{{{CO}}_{2}}) is the impulse-response function (IRF) for CO₂ that ranges from 1 at the time of the emission pulse (t = 0) to 0.41 after 100 years, and here it is set to a mean value of 0.52 over 100 years⁶⁰. The time horizon (TH) of our GWP calculations was fixed at 100 years following IPCC standards and previous studies^60,63,64.

Global mean surface air temperature change

We estimated the 100-year global mean surface temperature change for the estimated RF by adopting an equilibrium climate sensitivity (ECS), defined as the global mean surface air temperature increase that follows a doubling of pre-industrial atmospheric carbon dioxide (RF = 3.7 W/m²). Given a value of RF induced by a forcing agent, the temperature change is estimated as RF/3.7 × ECS. To consider the uncertainties of ECS, we adopted a mean value of 3 °C and a very likely (90% confidence interval) range of 2–5 °C following IPCC AR6⁵⁰. Without knowing the exact distribution shape of ECS and future albedo-change-induced RF, we created a log-normal distribution (Supplementary Note 4) to approximate their asymmetric distribution through numerical simulation. We then conducted Monte Carlo simulations that draw 5000 random samples from each distribution to jointly estimate the uncertainties of global mean surface air temperature changes. We report the mean and 90% interval ranges of the change in temperature.

Source: Ecology - nature.com