Enhanced surface urban heat islands due to divergent urban-rural greening trends

Satellite observations show that the surface urban heat island intensity (SUHII) has been increasing over the last two decades. This is often accompanied by an increased urban-rural contrast of vegetation greenness. However, the contribution of uneven vegetation trends in urban and rural areas to the trend of SUHII is unclear, due to the confounding effects of climate change and changes in man-made infrastructures and anthropogenic heat sources. Here we use a data-model fusion approach to quantify such contributions during the peak growing season. We show that the LAIdif (the urban-rural difference of leaf area index) is increasing (P< 0.05) in 189 of the selected 228 global megacities. The increasing trend of LAIdif from 2000 to 2019 accounts for about one quarter of the trend in satellite-derived SUHII, and the impact is particularly evident in places with rapid urbanization and rural cropland intensification. The marginal sensitivity of SUHII to LAIdif is the strongest in hot-humid megacities surrounded by croplands and in hot-dry megacities surrounded by mixed woody and herbaceous vegetation. Our study highlights the role of long-term vegetation trends in modulating the trends of urban-rural temperature differences.


Introduction
Today more than 55% of the world's population lives in urban areas and this number will reach 68% by 2050 (United Nations 2019). Urbanization is one of the most visible anthropogenic land use/land cover changes (Weng et al 2004, Peng et al 2012. In addition to producing anthropogenic heat fluxes, urbanization alters many other aspects of the surface energy budget, resulting in higher surface and near-surface temperatures in cities than their rural surroundings, a phenomenon known as the urban heat island (UHI) effect (Zhao et al 2014, Oke et al 2017, Manoli et al 2019, Li et al 2019a. With global warming (Hansen et al 2010), the combined effects of UHIs and heatwaves will pose significant threats to public health (Kovats and Hajat 2008, Li and Bou-Zeid 2013, Mora et al 2017, Liao et al 2018a and increase the global cooling power consumption (Li et al 2019b). Previous studies suggest that the spatial variations of the surface UHI intensity (SUHII), which is defined as the difference of land surface temperature (LST) between urban and rural areas, are strongly controlled by background climates and mainly follow the precipitation gradient, but the spatial variations of significantly correlated with the SUHII (Weng et al 2004, Oleson 2012, Peng et al 2012, Chakraborty and Lee 2019. Owing to land-use management and climate change, trends in vegetation greenness over urban and rural areas are often different, which can either enhance or reduce the urban-rural vegetation contrast, thereby altering the SUHII (Oleson 2012. Satellite observations show widespread enhancements in leaf area index (LAI) surrounding built-up areas . Meanwhile, many urban areas become comparably less green, or exhibit declines in vegetation greenness, and are expected to become more impervious (Seto et al 2012, Schneider et al 2015. Yao et al suggest that the global rural greening trend correlates with the increasing trend of satellite-derived SUHII (2019). However, it is unclear the extent to which the SUHII trend is directly caused by the biophysical effects of uneven vegetation changes between urban and rural areas (Chakraborty and Lee 2019), which is confounded by the simultaneous changes in climatic conditions, anthropogenic heat fluxes, and man-made infrastructures (Kalnay and Cai 2003, Zhao et al 2014, Li et al 2019a.
Here we aim to quantify the contribution of vegetation's long-term trend to the SUHII trend (trends are henceforth indicated by ∆). We study LST instead of near-surface air temperatures because LST is directly constrained by the surface energy balance equation (Chen et al 2020a) and its apparent change is observable from satellites. It has been reported that the SUHII is several folds stronger than the intensity of near-surface air UHI, making it easier to detect the vegetation biophysical effects (Venter et al 2021). We use a process-based attribution framework (the two-resistance mechanism method, TRM) to isolate the effect of LAI dif trend on the SUHII trend (Rigden and Li 2017, Liao et al 2018b, Li et al 2019a, Moon et al 2020, Chen et al 2020a, Wang and Li 2021) (see section 2 and supplementary information (available online at stacks.iop.org/ ERL/16/124071/mmedia)). This framework is based on the surface energy balance equation in which changes in vegetation greenness primarily alter the albedo for shortwave radiation, the aerodynamic resistance for scalar transfers, the surface resistance for water vapor transfer, the ground heat flux, and the emissivity (Chen et al 2020a). Acknowledging that the mechanism causing SUHII differ between daytime and nighttime (Peng et al 2012, Yang et al 2021 as a starting point we analyze the impacts of vegetation dynamics on the trends of daily average SUHII. We focus on the peak growing seasons, the three consecutive months with the largest mean LAI, during which the vegetation biophysical effects on LST are the strongest (Imhoff et al 2010, Peng et al 2012.

Land cover and definition of urban and rural patches
Land cover information is provided by the Collection 6 Moderate Resolution Imaging Spectroradiometer (MODIS) land cover product (MCD12C1, yearly, 0.05 • × 0.05 • ) (Friedl and Sulla-Menashe 2015). We use the International Geosphere-Biosphere Programme (IGBP) classification layer to define contiguous urban and rural areas (figure S1). First, a pixel is termed urban if it is classified as 'urban and builtup lands' by IGBP in any year between 2001 to 2018. This allows us to capture a maximum possible urban extent in the presence of urbanization, and is equivalent to assuming that a pixel can be transformed into an urban area and that the process is irreversible (Schneider 2012). A fixed urban extent also makes the trend of urban LAI unambiguous since allowing the urban extent to vary might cause an LAI trend even if the LAI of each pixel does not change. This implemented, a fixed urban extent effectively captures the underlying vegetation changes due to urban footprint changes (Yang et al 2019). Second, we label all adjacent urban pixels (i.e. eight-neighbours) as an urban patch. Each urban patch must have at least 16 pixels, which is about 493 km 2 at the equator. Third, in order to match the resolution of the reanalysis data (i.e. the enhanced land component of the European Centre for Medium-Range Weather Forecasts Reanalysis, or ERA5-Land), we aggregate the pixels within each urban patch from 0.05 • × 0.05 • to 0.1 • × 0.1 • resolution. During the aggregation process, a coarse-resolution pixel is marked as urban only if at least two of the four fine-resolution pixels are classified as urban (Schneider et al 2009). We repeat the second step to further merge the adjacent urban patches newly produced by the aggregation. Fourth, after identifying the urban patch, we define the rural patch as a surrounding 0.1 • -buffer zone (∼11 km), which allows the urban and rural areas to be approximately equal. Rural patches include only land pixels. Finally, we term a pair of urban and rural patches as a megacity. We analyze 228 megacities worldwide after excluding five megacities where the TRM framework fails due to negative aerodynamic or surface resistances (physically meaningless) (Liao et al 2018b). They are Cape Town (South Africa), Samarkand (Uzbekistan), Lahore (Pakistan), Almaty (Kazakhstan), and Dalian (China). Due to the spatial aggregation in our pre-processing of the land cover data, some megacities consist of multiple geometrically contiguous municipal cities, such as megacities of the Rhine-Ruhr (Germany), New York and Philadelphia (USA), the Greater Los Angeles (USA), the Pearl River Delta (China), the Yangtze River Delta (China), and the Tokyo Capital Region (Japan), etc. We note that the uncertainties in the MODIS land cover maps may propagate to the identification of urban and rural areas (Friedl et al 2002). Since our study focuses more on LST and LAI contrasts between urban and rural areas and their associated changes at the patch scale the accuracy of urban patch identification should be sufficient for this study.

MODIS LAI and LST
MODIS Collection 6 LAI products (MOD15A2H and MYD15A2H) are used in this study (Myneni et al 2015a(Myneni et al , 2015b, which are 500 m, eight day composites. Ground data have been used to verify the quality of the LAI products (Yan et al 2016a(Yan et al , 2016b, and these products have been widely used in previous work , Piao et al 2020. In this study, the LAI products are filtered with quality flags and composited to monthly averages using the number of days as the weight of each original eight day composite. We exclude pixels contaminated by clouds, aerosols, shadows, snow, and ice, and fill the gaps using the climatological mean monthly LAI during the study period (Samanta et al 2011. The monthly LAI is further aggregated to 0.1 • × 0.1 • and averaged to obtain the peak growing season mean. The peak growing season is defined as the three consecutive months with the highest mean LAI for each three month moving window (figure S2(a)). Then we calculate the average LAI for urban and rural patches (figure S3), as well as their difference (LAI dif ), for each megacity. Finally, LAI trends (∆LAI) are estimated using the Mann-Kendall test.
The MODIS LST products (i.e. MOD11C3 and MYD11C3) are used to compute SUHII for each megacity (0.05 • × 0.05 • , monthly composite, 2000-2019) (Wan et al 2015a(Wan et al , 2015b. We average the daytime and the nighttime LSTs, and use the average as a proxy for daily mean LST. According to the product user guide, the following quality filtering is applied: 'Mandatory QA flags' , 'Emissivity Error flag' and 'LST Error flag' have to be '00' or '01'; 'Data quality flag' has to be '0' . All MODIS LST data are converted to the same spatial and temporal resolution as the LAI data. We note that the MODIS SUHII is not used to calculate the impact of vegetation trend on the SUHII trend, which is diagnosed by the attribution method in section 2.5. The trend in MODIS SUHII is a result of vegetation biophysical effects, human-induced effects, and large-scale atmospheric changes (Zhou et al 2019).

ERA5 reanalysis
We use monthly averaged variables from ERA5-Land reanalysis (0.1 • × 0.1 • , diel average, 2000-2019) as inputs for the attribution framework described in section 2.5, including albedo (α), surface solar radiation downwards (S in ), surface thermal radiation downloads (L in ), surface net thermal radiation (L net ), surface latent heat flux (LE), surface sensible heat flux (H), and surface pressure (P s ) (Muñoz-Sabater 2019). Further, aridity index is calculated as the ratio of annual precipitation (Pr) to potential evapotranspiration (PET). ERA5-Land is a replay of the land component of the ERA5 reanalysis with an improved spatial resolution, including an elevation correction for the thermodynamic near-surface state .
In addition, we use monthly averaged air temperature and specific humidity at the lowest pressure level from ERA5 reanalysis (0.25 • × 0.25 • , diel average, 2000-2019) (Hersbach et al 2019). Depending on the surface pressure, the height of these atmospheric variables is about 40-100 m. This height is assumed to be above the surface roughness layer (Lee et al 2011), although in reality this assumption may be violated over some forest or urban areas (Oke et al 2017, Novick and Katul 2020, which may undermine the validity of Monin-Obukhov similarity theory in calculating the turbulent fluxes.

Population and Köppen-Geiger climate zones
The Gridded Population of the World, Version 4 (GPWv4): Population Count, Revision 11 at the resolution of 30 arc-second for two single years (2000 and 2020) are used. This dataset is produced by the Center for International Earth Science Information Network at Columbia University (CIESIN 2018). We also use the global map of Köppen-Geiger classification at 1 km resolution (Beck et al 2018). In this study, the global climates are divided into four main groups (figure S2(c)): A (tropical), B (arid/dry), C (temperate), and D (continental). None of the 228 megacities are located in the E (polar) category.

Diagnosing the impact of vegetation trend on the trend of SUHII
Changes in the SUHII are due to the combined effects of vegetation dynamics, changes in anthropogenic heat sources and man-made infrastructures, and large-scale climate change. We quantify the impact of uneven urban and rural vegetation trends on the trend of SUHII based on the TRM method. The TRM method was first developed to attribute the effect of land cover changes on LST using paired flux measurements or climate model simulations (Rigden and Li 2017, Liao et al 2018b, Li et al 2019a, Moon et al 2020, Wang and Li 2021. In this study, we attribute the vegetation-induced LST change (i.e. the urban-rural LST difference that is solely caused by LAI difference between urban and rural areas, denoted as SUHII LAI ) to five biophysical factors, namely, albedo (α), aerodynamic resistance (r a ), surface resistance (r s ), emissivity (ε), and the ground heat flux (which is computed as the residual of energy balance and denoted as REB) (Lawrence et al 2018, Chen et al 2020a), as follows: where ∂Ts ∂LAI is the sensitivity of LST to LAI over the reference (i.e. rural areas) (supplementary information), LAI u and LAI r are the LAIs in urban and rural areas, respectively, LAI dif is the urban-rural LAI difference. It is stressed that SUHII LAI is different from the satellite-observed SUHII that is affected by a variety of other influences as discussed earlier. Since the goal is to quantify the impact of vegetation trends on the trend of SUHII, we replace the LAI dif by the trend in LAI dif (i.e. ∆LAI dif ): where ∆SUHII LAI is the trend in the SUHII induced by the trend in LAI dif and the sensitivity is denoted as ∂SUHII ∂LAI dif .

Background LAI dif and SUHII
Overall LAI is higher in rural areas than in urban areas (i.e. LAI dif < 0, figure 1(a)). We divide the world into 11 subcontinental regions ( figure S2(b)). Megacities in eastern South America (ESA), eastern North America (ENA), and Oceania (OC) have relatively large LAI dif values (<−1 m 2 m -2 , table 1), but the largest LAI dif value is detected in the equatorial megacity of Singapore and Johor Bahru (-3 m -2 m -2 ). The LAI dif values are also large in coastal regions in eastern South America (ESA) and East Asia (EA) (figure 1(a)). In arid and semi-arid regions, the LAI dif values are close to zero or weakly positive (table 1, figure 1(a)). For comparison with LAI dif , the daily mean SUHII from Terra MODIS is shown ( figure 1(b)). The daily mean SUHII from Terra (overpass time: 10:30 a.m./p.m. local time) agrees with that of Aqua (overpass time: 1:30 p.m./a.m.), with Aqua SUHII having slightly larger magnitude (slope = 1.05, figure 1(d)). Considering the longer time series and less cloud cover in the morning, we use the daily mean SUHII from Terra for the following analysis. The SUHII is the strongest in North America (i.e. WNA and ENA) (table 1) and is also large in many coastal regions of ESA, EA, southeast Asia, and Europe ( figure 1(b)). Notably, we find urban cool islands (i.e. SUHII < 0 K) in arid places such as Saharan Africa and West Asia ( figure 1(b)). Statistically SUHII is negatively correlated with LAI dif (slope = −0.92, p < 0.001, figure 1(c)) and LAI dif explains 33% of the SUHII variance ( figure 1(c)).

Trends in LAI dif
Of the 228 megacities, 189 show a decreasing trend of LAI dif (becomes more negative) during the peak growing season, 90 of which are statistically significant (p < 0.05, figure 2(a)). In contrast, only 39 megacities show an increasing trend in LAI dif , of which only two are statistically significant (p < 0.05, figure 2(a)). The absolute and relative values of ∆LAI dif have similar spatial patterns (figures 2(a) and S4), which are nonetheless different from the spatial pattern of the background LAI dif (figure 1(a)). The strongest declines in LAI dif , which exceed −0.25 m -2 m -2 decade -1 (or 25% decade -1 ), occur in East Asia (EA, table 1) where strong greening trends are observed in rural China likely due to land-use management, along with declines in urban greenness due to rapid urbanization (Wang et al 2012. LAI dif is also declining rapidly in growing megacities such as the Gulf of Guinea in Africa, Southeast Asia, Sydney in Oceania, and Calgary and Edmonton in North America (figure 2(a)).
A decreasing trend of LAI dif may imply one of the following three scenarios: (a) rural greening is faster than urban greening, (b) rural browning is slower than urban browning, or (c) rural areas are greening while urban areas are browning. Our results indicate that 133 of the 228 megacities are greening simultaneously in both urban and rural areas but the rural areas are greening faster (uG-rG), suggesting a slowdown in terms of urbanization. Such megacities are concentrated in ENA, Europe, and extend into coastal regions in EA ( figure 2(b)). In addition, 52 megacities show urban browning and rural greening (i.e. uB-rG), reflecting a moderate rate of urbanization. These are distributed worldwide but are almost absent in Europe. Together these two categories (i.e. with negative LAI dif trends) include 185 megacities. Of the remaining 43 megacities, 37 show both urban and rural browning (uB-rB), which imply rapid urbanization, and six show urban greening and rural browning (uG-rB).

Trends in SUHII
We quantify the SUHII trend caused by the LAI dif trend (i.e. ∆SUHII LAI ) using the TRM  method ( figure 2(c)). Due to the strongest LAI dif trend, EA shows the largest ∆SUHII LAI (mean = 0.069 K decade -1 ), followed by SSEA (South and Southeast Asia, mean = 0.06 K decade -1 ), and SSA (sub-Saharan Africa, mean = 0.037 K decade -1 ) (table 1). These regions have experienced rapid urbanization and population growth over the past 20 years (table 1). Of the top 20 megacities with the largest positive ∆SUHII LAI globally (mean = 0.14 K decade -1 ), 15 are in China. For individual megacities, Hefei (China) and Chennai (India) rank 1st and 2nd in terms of ∆SUHII LAI (0.29 and 0.28 K decade -1 , respectively) (table S1) and both of them show a trend (∆SUHII LAI ) exceeding 20% decade -1 (table S1). This is due to rural (cropland) greening and urban browning (i.e. uB-rG) in both megacities, which is reflected in their higher ∆LAI dif (Hefei in EA: −0.31 m 2 m -2 decade -1 , Chennai in SSEA: −0.11 m 2 m -2 decade -1 ) compared to the average ∆LAI dif in their respective subcontinent regions (table 1). The opposite vegetation trends in urban and rural areas imply an increase in the urban temperature and a decrease in the rural temperature, which lead to an elevated SUHII. The spatial patterns of ∆SUHII LAI and the trend of SUHII (i.e. ∆SUHII) directly computed from Terra MODIS are in reasonable agreement, but the magnitude of ∆SUHII LAI is on average about one-fourth of ∆SUHII (figures 2(c) and (d)). The difference between ∆SUHII LAI and ∆SUHII highlights that uneven vegetation trends in urban and rural areas are not the sole contributor to the trend in SUHII (Oleson 2012. Evaluating the role of non-vegetation factors is out of the scope of this study, which requires, at the very least, data of anthropogenic heat fluxes and modeling the urban surface energy budget (Wang and Li 2021). But those are challenging and subject to large uncertainties (Allen et al 2011, Oleson 2012, Zhao et al 2014, Jin et al 2019, Zheng et al 2021. The contribution of anthropogenic heat flux to the SUHII could be 0.75-7.5 K according to a previous fine-resolution modeling study (Wang and Li 2021), which may explain the opposite signs between ∆SUHII LAI and ∆SUHII over the U.S. Great Lakes region (figures 2(c) and (d)). According to the U.S. Energy Information Administration, there is a strong declining trend in energy-related carbon dioxide emissions in the Great Lakes region from 2000 to 2018 (EIA 2021). As a result, the anthropogenic heat flux might be also declining and thus ∆SUHII is weaker compared to ∆SUHII LAI .

The sensitivity of SUHII to LAI dif
The sensitivity of SUHII to LAI dif is a marginal sensitivity under the current climatological states. It represents how much the SUHII is altered with a small change in LAI dif . This sensitivity has been shown to reasonably capture temperature changes caused by small perturbations, such as those caused by the two-decade vegetation trends shown above (Chen et al 2020a). It should be used with caution for understanding the SUHII which encodes large perturbations (Chen et al 2020b). The TRM method decomposes the sensitivity of SUHII to LAI dif (i.e. ∂SUHII ∂LAI dif ) into five contributors (figure S5). Surface resistance, which represents the ability to bring water to the surface, dominates ∂SUHII ∂LAI dif in 177 of the 228 megacities among different biophysical factors ( figure 3(b)). Thus, the spatial variability of LAI-induced SUHII changes at the global scale is primarily controlled by changes in the ability of land areas to evaporate water. This finding is consistent with the conclusion of previous studies (Bateni andEntekhabi 2012, Li et al 2019a), namely, vegetation changes mainly affect LST through the ET process. The ∂SUHII ∂LAI dif in the remaining 51 megacities are dominated by aerodynamic resistance. The megacities where ∂SUHII ∂LAI dif is dominated by aerodynamic resistance have significantly lower aridity index (Pr/PET) than those dominated by surface resistance (mean ± 1 standard deviation: 0.31 ± 0.27 vs. 0.51 ± 0.27, two sample t-test: p < 10 -5 ). Aerodynamic resistance represents the convective heat transfer efficiency affecting both sensible heat and latent heat fluxes. Many megacities where ∂SUHII ∂LAI dif is dominated by aerodynamic resistance are located in arid to semi-arid western United States (Zhao et al 2014), Saharan Africa, and West Asia (figures 3(b) and S2(c)) . In these megacities, changes in aerodynamic resistance caused by the vegetation trend are more likely to affect SUHII through sensible heat flux. However, for those megacities in the eastern United States where the climate is relatively humid (figures 3(b) and S2(c)), changes in aerodynamic resistance are more likely to affect SUHII through changing the latent heat flux (Chen et al 2020a). The contribution of LAI-induced albedo changes is on average about 6% to the ∂SUHII ∂LAI dif (figure S5), which depends on the sensitivity of lightextinction to LAI change within the vegetation canopy. As expected, contributions from other factors (i.e. ground heat flux and emissivity) are not significant (Chen et al 2020a).
We further group the ∂SUHII ∂LAI dif according to the Köppen-Geiger climate zones and biome types ( figure S2(c)). First, SUHII is most sensitive to LAI dif in megacities located in arid/dry areas surrounded by mixed woody and herbaceous biomes (W + H, figure 3(c)). In most of these areas, ∂SUHII ∂LAI dif is dominated by aerodynamic resistance ( figure 3(b)). Second, ∂SUHII ∂LAI dif increases from megacities surrounded by forests to megacities surrounded by croplands in a given climate regime (except for the arid/dry climate) (figure 3(c)). The higher sensitivity for megacities surrounded by croplands implies that for an equal change in LAI dif , intensification of rural croplands can increase SUHII more strongly. Incidentally, the LAI dif trends are large for cropland-surrounded megacities, especially in the continental and tropical climates (figure S6). As a result of the combined effects of higher sensitivities and larger LAI diff trends, 13 of the top 20 megacities with the largest positive ∆SUHII LAI are surrounded by croplands (table S1).
Why is ∂SUHII ∂LAI dif higher in megacities surrounded by croplands than those surrounded by forests (figure 3(c))? Due to the relatively low LAI in croplands on average (figure S6), the absorption of photosynthetically active radiation (PAR) is unsaturated. Therefore, compared to natural vegetation, an increase in crop LAI can absorb more PAR, which provides more available energy to trigger the opening of stomata and increases the evaporative cooling, leading to a larger ∂SUHII ∂LAI dif . Furthermore, croplands potentially have more water available for ET due to irrigation (D'Odorico et al 2020). These combined effects of canopy structure, radiation field, and water availability shape the higher aridity index (i.e. wetter) in croplands than other biomes across different climate zones (figure 3(c)). All of these are consistent with the finding that ∂SUHII ∂LAI dif is mainly controlled by surface resistance (figure 3(b)). For megacities surrounded by forests, the low ∂SUHII ∂LAI dif is due to their high background LAI in both rural and urban areas (figure S3). The high background LAI makes the biophysical factors insensitive to LAI due to saturation effects, especially for radiation and surface resistance.

Conclusion
In this study, we quantify the biophysical effects of vegetation dynamics on the contrast of LST between urban and rural areas over 228 megacities. Our results show that increases in LAI dif contribute to a quarter of the increasing trends in SUHII. The biophysical effects of vegetation trends on SUHII trends are determined by two major aspects. First, they depend on the magnitude of LAI dif trends, which is a joint effect of human land-use management and climate change such as warming and elevated atmospheric CO 2 , Piao et al 2020. Trends in LAI dif are most evident in megacities in Asia. Second, vegetation biophysical effects also depend on the sensitivity of SUHII to LAI dif, which is controlled by background climate and biome type. We find that this sensitivity is the strongest in hot-humid megacities surrounded by croplands and in hot-dry megacities surrounded by mixed woody and herbaceous vegetation. However, their dominant mechanisms are different. The former is mainly controlled by surface resistance which affects the latent heat flux, while the latter is mainly controlled by aerodynamic resistance which affects both sensible heat flux and latent heat flux depending on water availability. The high sensitivities suggest that continued cropland intensification and semi-arid woody vegetation encroachment in rural areas may enhance SUHII; while increasing urban vegetation in these megacities can reduce LAI dif and thus mitigate SUHII. The highest background LAI dif and SUHII are concentrated in megacities surrounded by forests in North and South America, but the sensitivity of SUHII to LAI dif is relatively low in these regions because of the high background LAI in urban and rural areas. Therefore, to reduce SUHII by the same amount as in megacities surrounded by croplands, megacities surrounded by forests would need to increase more urban greenness or invoke additional measures such as adopting white roofs and reducing anthropogenic heat emissions.

Data availability statement
The data that support the findings of this study are openly available.

Acknowledgments
C C and T F K are supported by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research, Reducing Uncertainties in Biogeochemical Interactions through Synthesis and Computation Scientific Focus Area. T F K acknowledges additional support from the U.S. Department of Energy Early Career Research Program Award #DE-SC0021023. D L is supported by the U.S. Department of Energy, Office of Science, as part of research in Multi-Sector Dynamics, Earth and Environmental System Modelling Program. We thank the National Center for Atmospheric Research (supported primarily by the U.S. National Science Foundation) for providing supercomputing resources. The authors declare no competing interests.

Author contributions
C C designed the study, analyzed the data, and drafted the manuscript. All authors interpreted the results and revised the manuscript.