Urbanization alters atmospheric dryness through land evapotranspiration

‘ Urban Dry (Wet) Islands ’ (UDI/UWI) represents microclimate change that impacts ecosystems and human well-being. However, causes of the UDI (UWI) phenomena are not fully understood due to the lack of empirical data. Here, we quantify UDI (UWI) using global observations of atmospheric humidity, evapotranspiration (ET), and land surface characteristics across 25 large urban agglomerations. We show that UDI (17) and UWI (8) are closely tied to local ET, global warming, and ‘ Urban Heat Islands ’ through intertwined linkages with water and energy balances. UDI is most pronounced in humid vegetated regions where mean urban-rural annual ET differences are as high as 215 mm, whereas UWI is found in arid regions or in climates with dry summers. We conclude that ET can be used as a single variable to explain emerging urban environmental changes. Our study supports a concerted strategy of restoring nature ’ s ET power as effective ‘ Nature-based Solutions ’ to mitigate the negative environmental effects of urbanization. npj Climate and Atmospheric Science (2023) 6:149 ; https://doi.org/10.1038/s41612-023-00479-z


INTRODUCTION
More than half of the world's population lives in urban areas and this number is projected to rise by 68% by 2050 1 .Compared with 2012, urban land uses are projected to double by 2030 globally, mostly in developing countries 2 .Permanent land use change through urbanization modifies the radiative, physiological, and aerodynamic properties of terrestrial ecosystems 3,4 .Urbanization is known to significantly elevate surface temperature or air temperature (Urban Heat Island, UHI) [5][6][7][8][9][10] , reduce latent heat (i.e., evapotranspiration) 11 , and increase storm water runoff 12,13 in general (Fig. 1).
Several factors have prompted recent interest in understanding the change in atmospheric dryness as concurrent with UHI.UDI (or UWI) has important implications for understanding cloud formation 23 , rainfall intensity 26 , climate change 26 , human bioclimates 27 and thermal comfort 27 , and urban design and planning 17 .In addition, there is evidence that wildland fires at the 'urban-forest interface' are linked to the drying trend of urban atmosphere 28,29 .Even a small change in atmospheric moisture availability may still be sufficient to alter latent heating, atmospheric circulation and moisture transport, and precipitation patterns [30][31][32] .Changes in atmospheric VPD that is affected by air temperature, is a key indicator of atmospheric moisture.VPD is increasingly recognized as an essential variable for assessing atmospheric water stress, heatwaves, plant photosynthesis and ecosystem productivity, and carbon sequestration in urban ecosystems [33][34][35][36] .
Long-term mean atmospheric humidity has been used for quantifying regional evapotranspiration (ET) 37 which is known to be one of three important sources of atmospheric moisture (i.e., antecedent atmospheric water vapor, lateral convective water vapor source, and local ET) 38 .The recent global 'greening up' due to reforestation and other factors (i.e., climate change, CO 2 fertilization, nitrogen deposition and recovery from natural disturbances) 39 has increased ET that impacted water resource availability 40 and precipitation patterns 31,41 .However, loss of vegetation or 'browning' was common in urban agglomerations worldwide, notably in southeastern China and the U.S.A. 10,12,13 .Although vegetation is well known to play an important role in the water and energy balances 41,42 , the significance of vegetation ET in linking land ecohydrological processes with urban atmospheric moisture or climate feedbacks is rarely considered in urban greening designs to mitigate extreme climatic events (e.g., heatwaves and storms) and hydrology (e.g., storm runoff and flooding) (Fig. 1).Recent extreme heatwaves and precipitation events around the globe exceeded prediction by most climate models and prompted renewed interest in land-atmospheric interactions 43,44 .Existing climate models may not adequately describe land-atmospheric interactions 43,44 , and therefore there is a large uncertainty in projecting future heat and precipitation patterns.In addition, traditional remote sensing products (i.e., MODIS) for vegetation and ET often exclude built-up areas.The role of ET in connecting urban climate, ecosystem health, and extreme hydrology 11,14,24,45 is increasingly recognized.
However, most previous studies on the impact of urbanization on local climate 8 or hydrology 46 have focused on the role of urban impervious surface 47 with little emphasis on vegetation ET, a key process that regulates latent and sensible heat fluxes as well as ecohydrology 6,8,[14][15][16] .For example, existing literature generally attributes UDI to global warming and UHI, both of which cause elevated air temperature and reduced relative humidity 48 .Although the cooling effects of ET on air temperature are well documented 5,49 , limited knowledge is available about the contribution of ET to explain changes in atmospheric humidity in an urban environment.There is no consensus on how land cover change, ecohydrology, air temperature, and atmospheric humidity/ moisture availability interact during urbanization 23,44,[50][51][52] .
Advanced knowledge of the close connections among vegetation greening (or browning), ecohydrology, and microclimate has practical applications for mitigating urban environmental change 53 .Over the past few years, a series of innovative strategies or practices have been explored to address environmental change challenges in cities. Examples include green infrastructure in Europe 54 , low-impact development in the United States, watersensitive urban design in Australia, natural infrastructure in Peru, Nature-based Solutions (NbS) in Canada 55 , and 'Sponge Cities' in China 56 However, there is a significant knowledge gap about the effectiveness of these measures, resulting in low confidence among stakeholders and political decision makers to act 53,57 .A theoretical framework is needed to understand the tradeoffs of all the approaches in order to best leverage resources and finances for implementing these strategies 53 .Ecohydrological studies on the role of ET in surface energy and water balances are needed for devising sustainable urban green infrastructure 58 .
This study assembled a large, spatially, and climatically diverse data set from around the world to quantify UDI (UWI) as defined by various atmospheric humidity indicators at multi-spatiotemporal scales.We expanded our previous studies [14][15][16] with the largest dataset with flux measurements in urban areas and explored the casual effects between air dryness and ecohydrology.The overarching goal of this study is to provide more evidence of UDI (UWI), and most importantly, to mechanistically explain how atmospheric dryness is controlled by ET, a critical process that has been previously overlooked in understanding land-atmospheric interaction and ecohydrology in an urban environment.Such knowledge is important to design any 'Nature-based Solutions' strategies as sustainable means to mitigate the negative effects of urbanization on both ecosystems and human well-being globally.

RESULTS AND DISCUSSION
The conceptual framework Here we present an integrated conceptual model as a framework for understanding the physical processes induced by land use change that controls UHI and UDI (UWI) processes (Fig. 1).In general, the land-atmosphere interactions and exchanges of water and heat mainly depend on two factors: the vertical exchange rate through turbulence and differences in humidity or temperature between the ground surface and the atmosphere above.Previous studies [6][7][8][14][15][16] show that converting vegetated land to 'dry' or 'brown' urban impervious surfaces alters the land-atmosphere exchanges of water and heat, and thus the hydrological and meteorological processes (Fig. 1). Proceses and variables associated with the water and energy exchange between vegetation and the atmosphere and regional climate feedback include soil heat capacity, vegetation leaf area, surface albedo, aerodynamic roughness, and net radiation.All these variables control latent heat (i.e., ET) and sensible heat transfers (Fig. 1).
Fig. 1 A conceptual model illustrating the physical urbanization processes that contribute to 'Urban Heat Island' and 'Urban Dry (Wet) Island' (rural on the left and urban on the right).a Focuses on surface energy balance, R a is radiation, R e is reflection, H is the sensible heat flux, LE is the latent heat flux, Q s is the heat storage, Q A is anthropogenic heat, and A H is advection heat from the countryside.EZ is Entrainment Zone, light blue and light red arrows indicate urban-rural breeze circulation.b Focuses on the surface water balance, Pre is precipitation, T is transpiration, E soil is soil evaporation, E wat is water evaporation, E int is interception evaporation, A Q is advection moisture from the countryside, E A is anthropogenic vapor injection due to combustion and etc., and Q is runoff.In (c), the black or red lines represent before or after urbanization, respectively.r a is aerodynamic resistance, r s is surface resistance, T s is surface temperature, T a is near-surface air temperature, q Ã s is surface saturated specific humidity that is only a function of T s , q s is surface specific humidity, and q a is atmospheric specific humidity.The BLH is the air layer from the ground up to the entrainment zone, and the CLH is the air layer from the ground to canopy/roof level.The green and red colors for BLH and CLH represent rural and urban environments, respectively.The plus and minus signs in parentheses represent an increase or decrease, respectively.
We use this conceptual model to illustrate the likely influences of biophysical factors on the differences of urban-rural near surface atmospheric air temperature (T a ) and specific humidity (Δq a ).We hypothesize that UHI is primarily due to differences in evaporative cooling as opposed to aerodynamic resistance 7,8 and UHI processes are connected to UDI (UWI) through hydrothermal coupling.
In addition, we hypothesize that change in atmospheric humidity (Δq a ) is controlled by following major processes: changes in surface temperature (ΔT s ) (i.e., UHI), near surface air temperature (ΔT a ), radiative forcing (ΔR n ) (i.e., albedo), ΔG (i.e., diffusive heat flux into/out of the surface, including anthropogenic heat flux), aerodynamics (Δr a , e.g., wind speed), and surface evaporation and transpiration resistance (Δr s , e.g., vegetation canopy conductance, leaf area).However, the contribution of these terms to explain Δq a likely varies due to the differences in background climate and land surface characteristics.For example, Δq a in humid regions characterized by relatively high vegetation cover fractions in urban areas might be more responsive to Δr s or ΔR n than that in arid/desert regions characterized by low vegetation cover fractions in urban areas where Δq a is most sensitive to change in air temperature ΔT a (Fig. 1c).

Enhanced UDI (UWI) intensity globally
The assembled published datasets consist of 20 global urban eddy flux sites (Fig. 2a) 59 and 140 paired weather stations (urban vs. rural) across China 15 (Fig. 3a).We define that UDI occurs when near surface atmospheric humidity in urban areas is lower than that in surrounding non-urban areas (i.e., Δq < 0 or ΔVPD > 0), otherwise we call it 'Urban Wet Island' (UWI, Δq > 0 or ΔVPD < 0).We found that more than half of the 20 urban eddy flux sites showed UDI effects (Δq < 0) in growing season (Supplementary Table 1, Supplementary Fig. 1), which were synchronized with UHI effects (Supplementary Table 1).
UDI effects were more pronounced in tropical and temperate regions than cold and arid ones (Fig. 2b).However, exceptions existed and UWI (Δq > 0) could occur in temperate climates such as CA-Sunset (California, USA) and GR-HECKOR (Heraklion, Greece).These two sites were dominated by a Mediterranean climate, typically having a 'dry summer' but exhibiting UWI effects.This dry summer UWI effect is similar to what is found in areas with an arid climate such as US-WestPhoenix in the earlier stage (USA) (Supplementary Fig. 1, Supplementary Table 1).In cold climates with no dry season or with dry winter only, such as US-Minneapolis (US) and KR-Ochang (Korea), characterized by relatively high vegetation cover fractions in urban areas, the urban q was generally greater than non-urban q in the summer (Supplementary Fig. 1).However, urban q in FI-Torni (Finland, Cold region), characterized by low vegetation cover fractions in urban areas, was generally lower than non-urban q in the summer (Supplementary Fig. 1).It appears that vegetation cover and  1) with four Köppen-Geiger climate regions, i.e., Temperate (with wet summer, TE), Cold (with wet summer, CD), Arid (or temperate with dry summer, AR), and Tropical (with rainforest, TR), respectively 85 .b Average annual Δq (1990-2020) and monthly Δq (2001-2020) between all 25 urban sites and non-urban areas in four climatic regions, respectively.Δq for other 20 global sites is difference in q between 20 urban sites (These data can be accessed from https://doi.org/10.5281/zenodo.7104984) 59and surrounding non-urban pixel (data from ERA5 as non-urban areas), respectively.Δq for 5 Chinese sites represent observed differences in average q between urban and rural sites in five urban agglomerations.The shading areas illustrate the mean and standard deviation (s.d.) of Δq.
background climate features play an important role in determining Δq.
We conducted a more detailed analysis using data from the 140 paired sites in China that complements the 20 global flux sites.In addition to atmospheric q, vapor pressure deficit (VPD) was used to identify diurnal, seasonal, and annual variations of air dryness.We found a similar pronounced UDI pattern (Δq < 0 or ΔVPD > 0) in temperate humid and sub-humid regions (Fig. 3a, b) and this pattern was amplified during daytime and in growing season (Supplementary Fig. 2, Supplementary Fig. 3).The nighttime UWI (Δq > 0, ΔVPD < 0) effects showed opposite trends (Supplementary Fig. 3).Both Δq and ΔT a significantly contributed to ΔVPD during the entire study period (Fig. 3a).However, the rise of ΔT a alone was not sufficient to explain the observed increase in ΔVPD for the  for all 20 urban sites worldwide.a Significant correlations between LE versus q for Δq < 0 and Δq ≥ 0, respectively.b Correlation coefficient of LE versus q for individual sites.All the 20 sites passed Two-tailed T-test at the significance level of 0.001.The names and locations of the 20 sites are in Fig. 1a and Supplementary Table 1.Δq was calculated using data from 20 urban sites published in https://doi.org/10.5281/zenodo.7104984 59, and surrounding pixel data from ERA5 as non-urban areas, respectively.LE and q data were monthly data (n = 621) from 20 urban flux sites (2003-2020).
post urbanization period (2001-2018) (Fig. 3b).Attribution analysis found that the overall contribution of Δq to ΔVPD was much greater during the later period 2001-2018 than that during the first period of 1980-2000, especially in daytime in humid and subhumid regions (Fig. 3a).The variable relationships between ΔT a , Δq, and ΔVPD over time offered additional evidence that stage of urbanization affected air dryness characteristics.

Urbanization affects atmospheric dryness through ET
Both UDI and UWI were found in each of the 20 global urban sites studied (Fig. 2a, b).Regardless of UDI or UWI, monthly Δq fluctuates in the growing season (Supplementary Fig. 1).UDI effects were common in cities in tropical and temperate regions with wet and warm summers (Fig. 2b).This corresponded to lower ET rates in urban cores that had a smaller fraction of vegetative cover compared to their rural counterparts in the humid regions (Fig. 2a, Supplementary Fig. 4).Lower ET rates in urban areas have been widely reported in literature for temperate and subtropical regions in the US 11,46 and China 14,60 The exceptions found in CA-Sunset, a UWI case (Supplementary Fig. 1) can be explained by the lower ET in rural areas during dry summerlawn irrigation might be used in urban areas.In wet and cold climates regions, however, ET is often rather low due to energy limitations in winters.In the warm summer growing seasons, ET rates from well-vegetated sites of US-Minneapolis1 and US-Minneapolis2 in northern U.S., KR-Ochang in Korea, and FI-Kumpula in Finland (Supplementary Fig. 1) were generally greater than rural areas due to relatively higher temperature and thus higher VPD in urban areas.Similarly, USB (Baltimore, US) and NLA (Amsterdam, Netherlands) showed UWI effects due to high tree area fraction (54%) and water area fraction (17%), respectively (Supplementary Fig. 1, Supplementary Table 2) 59 .PLL (Lipowa, Poland) and PLN (Narutowicza, Poland) had similar climate and environment, but Δq was opposite (UWI in PLL, UDI in PLN) possibly due to their different aerodynamic roughness length for momentum (0.94 m at PLL and 1.99 m at PLN), respectively, calculated by the Macdonald morphometric method 59 .However, compared with FI-Kumpula, FI-Torni in Finland had a lower urban vegetation cover, urban q was generally lower than non-urban q in the summer (Fig. 2, Supplementary Fig. 1) corresponding to less ET in urban areas.The correlations among annual average Δq and ΔT a , Bowen ratio, and biophysical drivers of impervious surface fraction, vegetation fraction, displacement height, and roughness length for momentum (Supplementary Table 2, Supplementary Fig. 5, Supplementary Fig. 6) indicate that differences in surface evapotranspiration resistance, as well as surface aerodynamic properties such as roughness and displacement height between urban and rural areas all contribute to observed UDI and UHI intensity.
Further analysis of the eddy covariance flux measurements revealed that UWI (Δ q > 0) or UDI (Δ q < 0), latent heat flux (LE) or Evapotranspiration (ET) has significant positive logarithmic correlations with atmospheric specific humidity (q) (P < 0.001) (Fig. 4a, b) at a monthly scale for most of the 20 urban sites.However, the correlation was stronger for UWI (Δ q > 0) than that for UDI (Δ q < 0) (Fig. 4a), and in those cities with higher vegetation fraction (p < 0.001).We also found LE was significantly positively correlated with vegetation fraction, but negatively with impervious surface fraction (Supplementary Fig. 7a, c) (P < 0.001) at urban sites.On the contrary, sensible heat (H) was significantly negatively correlated with vegetation fraction, but positively with impervious surface fraction (P < 0.001) (Supplementary Fig. 7b, d).The linkage of LE (or ET), vegetation, and specific humidity indicated that the decline of vegetation due to urbanization impacted atmospheric specific humidity through ET.
Using independent remote sensing-based ET data, we found that 19 out of the 20 urban eddy flux sites had much lower ET than surrounding non-urban areas (Fig. 5).The mean annual ET difference is about 215 mm.There were significant positive correlations between ΔET and Δq as well as significant negative correlations between ΔET and ΔT a in those cities with lower urban ET and synchronous occurrences of UDI (Δq < 0) and UHI (ΔT a > 0) (solid circles in Fig. 5, Supplementary Table 1).These findings signified the close connections among ET, UDI, and UHI.On the other hand, the observation based on LE and q at the 20 eddy flux urban sites showed that UDIs (i.e., Δq < 0) were generally found in those urban sites with significant positive correlation of LE vs. q, and with higher vegetation cover losses in urban areas (Supplementary Fig. 8).Therefore, both remote sensing data and in situ flux data confirmed the likely role of ET (i.e., LE) in affecting UDI intensity (Figs.4a, b, 5).Urban sites with UWI also showed significant positive correlation between LE and q, due to various factors (e.g., irrigation or artificial greening), being consistent with the hypothesis that ET controls q.
Similar to the q-ET patterns found across the 20 global urban eddy flux sites, mean annual ET rates over urban areas were generally lower than the surrounding non-urban areas in the temperate southeastern China where UDI was most pronounced (Figs. 2, 3).Under such climates, the decreased trend of annual ET in urban areas followed a similar pattern to leaf area index (LAI) (Supplementary Fig. 4).Similar to an arid site (US-WestPhoenix) in the global eddy flux dataset, cities in the arid northwest of China had higher q and ET than surrounding rural areas (i.e., UWI) (Figs. 2, 3a).Generally, urban ET is lower than non-urban ET in energy-limited regions, while urban ET is higher than non-urban ET in moisture-limited regions.In humid energy limited regions where soil water does not limit ET, biophysical factors such as leaf biomass that directly affect ET become more important 15 .In addition, extensive impervious surfaces also limit land evaporation 14 .In arid/desert regions in western China, evaporative energy is not limited due to the abundant incoming solar radiation and The size of circles represents magnitudes of correlation coefficients of average monthly ΔET vs Δq, ΔET vs ΔT a , respectively (see Supplementary Table 1).ΔT a (°C), Δq (g/kg), and ΔET (mm), were the differences in near-surface air temperature, specific humidity, and evapotranspiration between urban and surrounding non-urban areas.Δq and ΔT a were calculated using 20 urban sites dataset from https://doi.org/10.5281/zenodo.7104984 59, and surrounding pixel data from ERA5 as non-urban areas, respectively.ΔET pixel data are from remote sensing SSEBop ET products 81 .The information for the 20 urban flux sites is in Fig. 1a.Confidence levels are denoted on labels by *P < 0.001 and + P < 0.01.high VPD.Urbanization does not necessarily mean a decrease in LAI in arid regions where artificial greening using irrigated plants is common resulting in oasis 5,15,61 .
The temporal variations of ET and q during urbanization over China (Fig. 3) offered additional evidence about how the two variables interact.Annual ET rates generally decreased during 2011-2017, the accelerated urbanization period, compared with 2003-2010, the initial urbanization period, in most urban cores except the arid TSB in western China (Supplementary Fig. 4).The most significant decrease in ET occurred in eastern China.Similar to results from analysis performed over urban cities worldwide (Fig. 4), ET in cities in eastern China also positively correlates with q, which indicates that the reduction in ET as a result of urbanization reduced the content of water vapor in the atmosphere (Supplementary Fig. 4).The consistent spatiotemporal variation patterns of atmospheric humidity q, air temperature T a , ET and LAI (Fig. 3, Supplementary Fig. 4) in southeast China illustrate the close coupling of these climatic and ecohydrological variables.
Global data from both eddy covariance flux measurements and remote sensing products all reveal that urbanization impacts atmospheric dryness and air temperature through altering ET processes (Figs. 4, 5, Supplementary Fig. 7, Supplementary Fig. 4).Loss of vegetation in humid and sub-humid cities or irrigation and artificial planting in arid regions alters the relationship between T a , q, and ET (Supplementary Fig. 4).Cities could get hotter and drier than surrounding rural areas in most temperate or humid regions because of a decline in ET (Figs. 3, 4, Supplementary Fig. 4).The UDI pattern is amplified during daytime (Fig. 3a, Supplementary Fig. 2) and in growing season (Fig. 2, Supplementary Fig. 3) due to the large contrast in ET between urban vs. rural areas.

UDI (UWI) and UHI coupling through ET
The 'surface flux equilibrium' theory 37 suggests that the variables of near-surface air temperature T a and specific humidity q are sensitive to turbulent surface fluxes, sensible heat flux H and latent heat flux or ET, respectively.The theory implies that urbanization that generally leads to a decrease in ET and an increase in the ratio of sensible heat exchange to latent heat exchange (i.e., the Bowen ratio) would decrease air humidity and increase air temperature.
Our analysis on monthly urban surface fluxes in four climate zones at the 20 urban flux sites confirmed this hypothesis with exceptions in extreme conditions (Fig. 6).Lower Bowen ratio (B ≤ 1) was generally associated with higher LE (ET), leading to a wetter (increased q) near-surface atmosphere, while higher Bowen ratio (B > 1) was associated with lower LE, leading to lower q in urban areas, such as in temperate and cold climate (Fig. 6a, b).Exceptions were found in arid regions where LE is not associated with lower Bowen ratios (Fig. 6d).Another exception found was in regions under a tropical climate where LE was similar (Fig. 6c).Nevertheless, our data show that vegetation cover and background climate play an important role in determining local q through ET on a monthly time scale.
In addition to ET-q connections under different Bowen ratio conditions, relative humidity and VPD are directly tied to air temperature, thus UDI (UWI) and UHI are coupled.Temperature alone may not provide a complete explanation of the heat storage changes in the Earth system 38 for evaluating future global climate change.ET is significant to explain the integrity of urban physical process of hydrothermal coupling.The analysis by Li et al. 8 and this study confirmed the important role of ET in regulating both surface temperature and atmospheric dryness.The close linkages among ecohydrology (i.e., ET), UDI (UWI), and UHI are of great significance to modeling the entire hydrothermal coupling processes in urban ecosystems 24,62 .Therefore, the UHI and UDI phenomena should be collectively addressed in landscape design, urban adaptive strategies development, and global change assessment 14 .
Vegetation ET: key to understanding urban environmental change Evapotranspiration is a keystone climate variable that uniquely links the water cycle (evaporation), energy cycle (latent heat flux), and carbon cycle (transpiration-photosynthesis trade-off) 8,63,64 .Forested watersheds have lower stormflow than urban ones because the former have higher ET in the growing season reducing antecedent soil moisture and increasing water storage 46 .Urbanization exacerbates global warming and UHI effects on atmospheric drying through reducing ET and atmospheric water availability 15,16 .
The present and previous studies 16 show that UDI is more pronounced during the daytime in humid regions and during growing seasons, thus air dryness is closely linked to plant ET dynamics.On the other hand, ET in arid cities can be higher than in rural areas due to artificial urban greening and irrigation, especially during drought conditions 11 .Overall, urbanization is more likely to cause 'Oasis Wet Islands' (Δq > 0) effects in the arid zone (e.g., the TSB area in western China and California sites in the dry summer), while triggering 'Urban Dry Island' (Δq < 0) effects in the humid regions.Left is a lower Bowen ratio (B ≤ 1) and right is a higher Bowen ratio (B > 1).q is specific humidity near-surface atmosphere, T a is air temperature near-surface atmosphere.H is sensible heat flux and LE is latent heat flux.LE, H, q, and T a were measured at all 20 urban sites by the eddy covariance method 59 , and were standardized for comparison.Size of arrows represents relative magnitudes.
Recently, the scientific community has shown great interest in how land cover change affects or downwind precipitation, moisture availability, and water supply 31,41 .The contribution of local evapotranspiration, which could be altered by changes in land surface characteristics (e.g., deforestation or reforestation), to precipitation and moisture 'recycling' is therefore of considerable interest 30,65 .As one of the most important anthropogenic influences on land use/land cover change, urbanization has been found to increase mean and heavy precipitation over and/or downwind of cities with medium confidence 66 .However, there are limited studies about the magnitude the contribution of ET and surface condensation in urban areas 17 .In addition, quantification of the influence of ET on precipitation remains particularly uncertain 67 .Our process-based analysis on field data revealed that ET does affect the specific humidity at a monthly scale (Fig. 4a, b; Fig. 6), which contributes to the understanding of land surface (ET)-atmosphere (humidity) interactions, one of the least developed research areas in global climate change studies that primarily rely on mathematical modeling, rather than observations.
Our results contradict to certain modeling work that suggests that near-surface atmospheric humidity is largely controlled by the ocean and atmosphere, with the land surface playing a relatively passive role 68 .We argue that this discrepancy arises due to spatial and temporal scales concerned.At a long time and global spatial scales, the ocean and atmospheric circulations dominate air humidity 37,69,70 .However, at a monthly time scale, atmospheric specific humidity is closely associated with ET, a variable that is highly influenced by local land use and land covers.Our study used 'paired' method to detect the sole local effect of urbanization.This 'paired' method removed the effects of ocean and lateral convective heat and energy transfers, and thus provided a solid evaluation of the 'net' urbanization effects on local climate.However, other physical processes such as anthropogenic vapor sources due to combustion and artificial irrigation, convective flux of heat and moisture via entrainment, and advection moisture from surrounding areas that contribute to the UDI/UWI were not explicitly examined in this study.Other limitations may include site representations and observational variables.Although the 20 global urban sites are reasonably distributed across mean temperature and precipitation for global urban locations, data gaps remain 59 .
Currently, the eddy flux network and ET measurements for urban conditions are sparse 11,44 and remote sensing-based ET products are rarely validated for urban areas 63 .The overall dominance of evapotranspiration for determining variability of UDI intensity in our study is consistent with previous studies in temperate humid region YRD and PRD, China 14,15 and in temperate humid city Boston, U.S.A. 24 .Overall, our studies offer strong empirical evidence on the close land surface (ET)atmosphere (humidity) interactions.

Implications for urban natural-based solutions
One of the many Nature-based Solutions measures 71 for mitigating detrimental urbanization effects (i.e., UHI, heatwaves, or storm water runoff) is urban forestry: tree planting 72 or installing other green infrastructure (i.e., green roofing, wetland restoration) in urban watersheds 12,13,53,73 .This study quantitatively explains the fundamental linkages between atmospheric humidity, a major driver of urban ecosystem productivity, land surface processes, and ecohydrology (ET).Such a framework provides the basic science and guidance for designing optimum greener, cooler, and spongy cities.The attribution analysis results support our previous recommendations that mitigating the negative effects of UDI, UHI, storm water runoff should focus on restoring the ET capability of urban ecosystems as an effective means of 'Nature-based Solutions' (NbS) 15,58,60 .
The diverse synergies and trade-offs among ecosystem services associated with NbS have important implications for urban planning and management under different climates 53,74 .For example, urban greening may lead to elevated ET and humidity which may contribute to concerns of increasing heat index and human heat stress, thus may counteract the benefit gained from evaporative cooling in humid areas 74 .Because UHI is controlled by urban-rural differences in ET and convection efficiency, increasing green cover and albedo are more efficient in dry regions than in tropical cities to achieve cooling effects 75 .However, green roofs and parks are more effective in achieving cooling potential in wet regions than in arid regions 53 .The arid regions can have both benefits of cooling and flood retention by city greening aided with irrigation if possible.Concerns about sustainable vegetation coverage and the negative impact on local water resources make the NbS more changing in dry areas.UDI (or UWI) can also synergistically interact with UHI and increase urban ecosystem risks.For example, increased atmospheric VPD (UDI) and air temperature (UHI) have been identified as important threats to urban tree photosynthesis, growth, and health 76 , as well as potential cause of wildland fires at the 'urban-forest interface' 28,29 .
Using global ET data of in situ eddy covariance measurements and remote sensing products, we confirm that, like other emerging environmental changes such as UHI and elevated stormwater runoff in human-dominated ecosystems, atmospheric moisture has been substantially altered in many parts of the world due to urbanization.The magnitudes and directions of changes depend on background climate and changes in land surface characteristics such as vegetation cover leaf area, and timing of the assessment (i.e., daytime vs. night, growing season vs. dormant season).Therefore, as the main regulators of evaporative cooling and atmospheric moisture in cities, vegetated ecosystems including urban forests and lawns, and natural or man-made wetlands should be carefully protected, managed, and designed as part of the urban green infrastructures.The trade-offs or synergies of various NbS for coping with climate warming, UHI, storm runoff, and human welfare must be recognized.The efficiency of UHI and UDI mitigation strategies varies across geographic and climatic regions, thus, any NbS efforts aiming at regulating city environmental change should be put in the context of local ecohydrological-climatic conditions.

Global databases on urban atmospheric moisture
We defined UDI (or UWI) and UHI intensity as the difference between the atmospheric humidity (q or VPD) and air temperature in the urban near-surface and the corresponding height for the undisturbed rural areas, respectively.The urban near-surface was chosen because of its close association with urban land use change and its importance to public health 74 .We assembled a large dataset that consists of ET, meteorological variables (e.g., q, T a , VPD, LE, H), and biophysical parameters from 20 urban eddy flux sites across thirteen countries 59 , and 140 paired weather stations in five urban agglomerations across mainland China, representing a wide range of physiographical regions from the boreal to the tropics.
For the global dataset, the specific humidity q and air temperature T a at non-urban areas were obtained from ERA5 reanalysis climatic data at 0.25°spatial and hourly temporal resolutions from 1990 to 2020 77 .The meteorological factors, such as specific humidity q, air temperature T a , latent heat flux LE, sensible heat H, and biophysical factors, such as impervious surface fraction, vegetation fraction, roughness length for momentum, and displacement height at 20 urban sites were obtained from https://doi.org/10.5281/zenodo.7104984 59(See Supplementary Table 3, Supplementary Tables 4, 5) (More details see Supplementary Methods).
For the dataset across China, we assembled long (1980-2018) daily climate data from 580 weather stations that represent the core meteorological monitoring stations in the five urban agglomerations in China and cover a large climatic gradient (Supplementary Fig. 9).The 6-hourly (02:00, 08:00, 14:00, and 20:00 BJT) meteorological observation data from 580 standard weather stations were acquired from the China Meteorological Administration (CMA) (http://data.cma.cn/en).We used near-surface (a height comparable to be measured at 1.5 m above the ground at standard weather stations) data to calculate atmospheric humid indicators, i.e., specific humidity q, vapor pressure deficit VPD and estimate potential ET (PET) 78 .A daily mean temperature and other variables were calculated by averaging four 6-hour observations in a day.

Identifying urban-rural pairs
The global urban-rural pairs were constructed by matching the 20 urban flux sites with ERA5 reanalysis climatic data (1990-2020) 77 .The land surface scheme used to develop the ERA5 products does not include urbanization effects 59,77 , representing the rural baseline climate 74 .Therefore, we considered the 20 urban flux tower sites (gap filled from ERA5) to represent urban climate and ERA5 data with 20 matched grids (0.25°× 0.25°) surrounding 20 urban flux tower sites as paired rural sites.
For selecting urban-rural pairs across China, we assembled long term (1980-2018) daily climate data from 580 weather stations that represented the core meteorological monitoring stations in the five urban agglomerations and covered a large climatic gradient in China (Supplementary Fig. 9).Using the Paired Urban-Rural Classification method 6,79 , we selected 140 paired urban vs. rural meteorological stations across China.More details see Supporting Information.The impervious surface dataset 80 with high spatial details and overall accuracies of more than 90% (http://data.ess.tsinghua.edu.cn/gaia.html)was used to identify paired urban-rural stations in China (Supplementary Fig. 9) (More details see Supplementary Methods).

Quantifying differences in Annual ET (ΔET) between urban and non-urban areas
The remote sensing-based SSEBop ET dataset 81 has the advantages over other products in that it covers urban areas.Therefore, we used this ET product to compare ET rates between urban and surrounding non-urban areas and establish a causal linkage between UDI (UWI) and ET (see Supplementary Table 4).Monthly SSEBop ET estimates at 1 km spatial resolution (https://earlywarning.usgs.gov/fews/product/458) was obtained from the operational simplified surface energy balance (SSEBop) modelling approach 81 .The impervious surface dataset (http://data.ess.tsinghua.edu.cn/gaia.html) 80with high spatial details was used to identify the border of urban worldwide (using four years of 2005, 2010, 2015, and 2018) in Fig. 2a.Referring to the UHI studies 82 , a 10 km buffer was created around all the urban areas to represent 'surrounding non-urban' landscapes (Fig. 2a).Then over the identified areas, monthly ET in urban areas and non-urban areas were estimated respectively and global ΔET was calculated using SSEBop ET dataset (2003-2020).

Determining biophysical factors controlling UDI (or UWI)
Correlation analyses were made for key variables at both 20 urban sites worldwide and pixel level across China including the pairs of ET and LAI, ET and q, ET and T a , Δq and vegetation fraction, Δq and Surface roughness etc.
The Global Land Surface Satellite (GLASS) LAI data sets (version 3.0) generated from Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance data at an 8-day and 1 km resolution were used to assess vegetation change and its correlation with ET from 2001 to 2017 in China.The accuracy of estimated LAI dynamics was extensively validated for the study period 83 .
The ANUSPLIN model was applied to interpolate the meteorological grid data with a 1 km spatial resolution.We used the partial thin-plate splines with a linear dependency on elevation, which was generated from digital elevation model (DEM) with a 1 km resolution (www.geodata.cn).
The nonparametric Mann-Kendall test was performed to detect the parameter trends.Five types of trends, 'Significant increase', 'Increase', 'No change', 'Decrease', and 'Significant decrease' for each of the pixels across China, were determined according to their Z values.
Considering 2000 was the most obvious urbanization turning point in China due to their significantly higher ISAs 15 , we separated the entire dataset (1980-2017) into two periods, the 'pre-urbanization period' (1980-2000) and the 'post-urbanization period' (2001-2018), to detect changes in q over time.
Determining factors controlling ΔVPD between urban and non-urban areas A more detailed analysis was conducted for the Chinese dataset to quantify the relative contributions of urban-induced changes in near-surface air temperature (ΔT a ) versus specific humidity (Δq) on the VPD gradient (ΔVPD), we used the attribution method proposed by Huang et al. 15 to identify the attribution of ΔVPD induced by ΔT a and Δq.Based on the Clausius-Clapeyron equation, the variability of VPD is directly affected by the variabilities of e s and e a .To quantify the contributions from H 2 O and T a spatial variabilities to the VPD spatial variability, we separate the change of VPD (dVPD) into two parts: the contributions from water vapor partial pressure variability (dVPD q ) and those from T a variability (dVPD Ta ).Applying the chain rule to dVPD, we get dVPD ¼ de s À de a (1)   Where e a is the water partial pressure, and e s is the saturation vapor pressure calculated based on the equation from Murphy & Koop 84 .We assume that the magnitudes of de s and de a terms are relatively small, such that only the first terms of the Taylor expansion are considered in Eq. 2 and 3: de a % Δðe a Þ ¼ e a;1 À e a;2 We also separate the changes in e a term de a into the contributions from the specific humidity (q, kg kg −1 or g kg −1 ) and the variability of total pressure (p, hPa).Similarly, only the first terms of the Taylor expansion are considered for d q and d p , assuming that changes are relatively small.
Where ε is the ratio of molecular weight of water vapor to air, which can be calculated from molar mass of water vapor M v and molar mass of dry air M d (ε ¼ M v =M d ¼ 0:622).

DATA AVAILABILITY
The

Fig. 2
Fig. 2 Temporal and spatial distributions of the contrasts of mean atmospheric specific humidity (Δq) between urban and surrounding non-urban areas in 25 major urban agglomerations.Background ΔET map is created using high resolution SSEBop ET data (2003-2020) 81 scaled to the provincial administration level.Cool background blue colors indicate that ET in urban is lower than that in surrounding nonurban areas.a Mean Δq (g/kg) at 25 urban sites (circles) during 2001-2020.Red and blue circles indicate Δq < 0 (i.e., Urban Dry Island, UDI), otherwise Urban Wet Island, the size of circle indicates the intensity of Δq.The labels are site name abbreviations (for full names see Supplementary Table1) with four Köppen-Geiger climate regions, i.e., Temperate (with wet summer, TE), Cold (with wet summer, CD), Arid (or temperate with dry summer, AR), and Tropical (with rainforest, TR), respectively85 .b Average annual Δq (1990-2020) and monthly Δq (2001-2020) between all 25 urban sites and non-urban areas in four climatic regions, respectively.Δq for other 20 global sites is difference in q between 20 urban sites (These data can be accessed from https://doi.org/10.5281/zenodo.7104984)59and surrounding non-urban pixel (data from ERA5 as non-urban areas), respectively.Δq for 5 Chinese sites represent observed differences in average q between urban and rural sites in five urban agglomerations.The shading areas illustrate the mean and standard deviation (s.d.) of Δq.

Fig. 3
Fig. 3 Impervious surface area (ISA) (2017) distribution and change in urban-rural differences in near-surface air temperature (ΔT a , °C), atmospheric specific humidity (Δq, g/kg), and their contributions to differences in vapor pressure deficit (ΔVPD, hPa) across China.a The bars illustrate relative individual contributions of Δq and ΔT a to annual mean ΔVPD diurnally across five urban agglomerations during 1980-2000 (pre urbanization period, bars in the first and third columns) and 2001-2018 (post urbanization period, bars in the second and the fourth columns).Contribution rates (%) of ΔT a are shown in the upper panel and contribution rates of Δq in the lower panels.The attribution methods of ΔVPD see Methods.Error bars are ± 1 s.e.m Confidence levels are denoted by * P < 0.01 and ** P < 0.001.Five urban agglomerations including the Pearl River Delta (PRD), the Yangtze River Delta (YRD), Chengdu-Chongqing (CY), Beijing-Tianjin-Hebei (JJJ), and the Northern Tianshan cities (TSB).b Three lines illustrate annual changes in Δq, ΔT a , and ΔVPD in PRD and YRD regions with significant UDI effects.

Fig. 4
Fig.4Measured mean monthly latent heat flux (LE) or Evapotranspiration (ET) by the eddy covariance method and specific humidity (q) for all 20 urban sites worldwide.a Significant correlations between LE versus q for Δq < 0 and Δq ≥ 0, respectively.b Correlation coefficient of LE versus q for individual sites.All the 20 sites passed Two-tailed T-test at the significance level of 0.001.The names and locations of the 20 sites are in Fig.1aand Supplementary Table1.Δq was calculated using data from 20 urban sites published in https://doi.org/10.5281/zenodo.710498459, and surrounding pixel data from ERA5 as non-urban areas, respectively.LE and q data were monthly data (n = 621) from 20 urban flux sites (2003-2020).

Fig. 5 A
Fig. 5 A summary of the relationships among average monthly ΔET, ΔT a , and Δq, respectively, for 20 urban sites (2003-2020).The size of circles represents magnitudes of correlation coefficients of average monthly ΔET vs Δq, ΔET vs ΔT a , respectively (see Supplementary Table1).ΔT a (°C), Δq (g/kg), and ΔET (mm), were the differences in near-surface air temperature, specific humidity, and evapotranspiration between urban and surrounding non-urban areas.Δq and ΔT a were calculated using 20 urban sites dataset from https://doi.org/10.5281/zenodo.710498459, and surrounding pixel data from ERA5 as non-urban areas, respectively.ΔET pixel data are from remote sensing SSEBop ET products81 .The information for the 20 urban flux sites is in Fig.1a.Confidence levels are denoted on labels by *P < 0.001 and + P < 0.01.

Fig. 6
Fig. 6 Urban surface flux equilibrium on a monthly time scale in four climate zones to show coupling between energy balances and air temperature and humidity.a Temperate climate.b Cold climate.c Tropical climate.d Arid climate.Left is a lower Bowen ratio (B ≤ 1) and right is a higher Bowen ratio (B > 1).q is specific humidity near-surface atmosphere, T a is air temperature near-surface atmosphere.H is sensible heat flux and LE is latent heat flux.LE, H, q, and T a were measured at all 20 urban sites by the eddy covariance method59 , and were standardized for comparison.Size of arrows represents relative magnitudes.