Enhanced drought resistance of vegetation growth in cities due to urban heat, CO₂ domes and O₃ troughs

Terrestrial ecosystems assimilate ca. 30% of anthropogenic carbon dioxide (CO₂) emissions (Friedlingstein et al 2020) and have been a substantial carbon sink over the past six decades (Ballantyne et al 2012, Keenan and Williams 2018, Morecroft et al 2019). Increasing land carbon uptake has been partly attributed to enhanced vegetation growth in response to rising atmospheric CO₂ concentrations (Schimel et al 2015, Zhu et al 2016), known as the CO₂ fertilization effect. However, sustained global warming, coupled with extreme climate events such as droughts, have the potential to offset this CO₂ fertilization benefit, causing depressions in vegetation productivity and even a net release of CO₂ into the atmosphere (Lewis et al 2011, Choat et al 2018). Therefore, revealing the interactive and competing effects of altered temperature (ΔT) and CO₂ (ΔCO₂) on the response of vegetation growth and productivity to drought (D; hereafter denoted as ΔT | ΔCO₂ | D effects) is critical for understanding terrestrial carbon cycle dynamics and for informing reliable mitigation and adaption plans for climate change (Reichstein et al 2013, Brodribb et al 2020).

Examining the ΔT | ΔCO₂ | D effects on vegetation growth is generally challenging although there have been some manipulative experiments (Apgaua et al 2019, Birami et al 2020). With controlled experiments, there are uncertainties in translating experimental results into ecologically realistic predictions of how vegetation will respond to ΔT and ΔCO₂ under drought (Dusenge et al 2019, Brodribb et al 2020) at a large scale. In addition, the limited coverage of species and environmental conditions in manipulated experiments may result in inconclusive findings (Duan et al 2013, Birami et al 2020). For instance, ozone (O₃) exposure can induce reductions in vegetation growth (Gregg et al 2003, Ainsworth et al 2012) but is often overlooked in manipulative experiments. While large-scale free air CO₂ enrichment (FACE) and open top chambers experiments may allow exposure of vegetation to, besides rising CO₂, various stress factors such as drought and elevated temperature simultaneously (Ainsworth and Long 2021), such facilities are typically costly to maintain in the long term (Calfapietra et al 2010). It is also difficult in imposing more than two global change factors in these manipulative experiments. As such, there is a dearth of consistent long-term observations for advancing our knowledge of the interactive and combined effect of ΔT | ΔCO₂ | D effects on vegetation growth. For example, stomatal closure in response to drought can result in low leaf intercellular CO₂ and an enhanced sensitivity of photosynthetic rate to rising atmospheric CO₂, i.e. 'the low intercellular CO₂ effect', thus leading to relatively larger carbon uptake under drought and elevated atmospheric CO₂ (Kelly et al 2016). This beneficial effect was demonstrated using two Eucalyptus species in short-term experiments where leaf area index remained unchanged while long-term acclimation to drought counteracted this benefit (Kelly et al 2016, Jiang et al 2021).

Here, we approach the interactive effects of environmental drivers on vegetation growth and productivity along an urban-rural gradient (Calfapietra et al 2015) from 2001 to 2018 in the conterminous United States (CONUS). This urban plant physiology concept advocates the use of urban-rural gradients as a test bed and as a cost-efficient means for plant physiological studies since cities are experiencing conditions, such as temperature and CO₂ enhancements, decades ahead of projected change for natural ecosystems (Zhao et al 2016, Wang et al 2019). In addition, global climate change has exacerbated drought in both frequency and severity in urban and its neighboring rural regions (Güneralp et al 2015, Vicente-Serrano et al 2020), posing tremendous pressure on urban-rural ecosystems (Zhang et al 2019). With the long-term monitoring of drought conditions (Dai 2013, 2021), urban-rural contrasts in environmental conditions such as surface O₃ gradient (ΔO₃, typically higher O₃ concentrations in rural regions) can be exploited to reveal the response of vegetation growth and productivity to multiple altered environmental factors under drought. In this study, we analyze 75 urban-rural pairs of various sizes among three climatic zones (arid, temperate, and continental climate zones; figures 2, S1, and table S1 (available online at stacks.iop.org/ERL/16/124052/mmedia)) for revealing environmental drivers underlying ecosystem response to drought. The monthly enhanced vegetation index (EVI, 1 km) (Huete et al 2002) was used as a proxy for natural (non-crop) vegetation growth and productivity during the growing season, i.e. June, July, and August (JJA). EVI values were spatially averaged for urban and rural extents individually over pixels of plant functional types (PFTs) commonly identified in each urban-rural pair. We utilized a linear mixed-effects model to remove variation in EVI induced by climate variables such as solar radiation so that response of vegetation growth to drought and non-drought conditions could be inferred from EVI anomalies (i.e. observations—linear model predictions, denoted as ${\text{EV}}{{\text{I}}_{\text{a}}}$). Drought resistance of vegetation growth ($\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$) was computed as the difference in ${\text{EV}}{{\text{I}}_{\text{a}}}$ between drought and non-drought conditions for each urban-rural pair. The objectives of this study are to reveal the urban-rural differences in drought resistance of vegetation growth (hereafter denoted as ${\text{dEV}}{{\text{I}}_{\text{a}}}$) and to understand the underlying drivers (e.g. irrigation) of such differences across the CONUS, particularly the role of ΔT, ΔO₃, and ΔCO₂, in explaining the observed ${\text{dEV}}{{\text{I}}_{\text{a}}}{ }$at a large scale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1. Data for vegetation, weather, and drought

2.2. Drought response of vegetation growth and productivity

Prior to the quantification of drought response of vegetation growth and productivity, common PFTs were identified (further details related to identification of common PFTs can be found in supplementary section 2). Albeit the isolation of common PFTs within urban-rural pairs, EVI signals at 1 km still suffer background spectral contaminations from impervious surfaces or other non-vegetation pixels, resulting in smaller EVI values in urban extents compared to those in rural areas (Zhao et al 2016) with exceptions for urban-rural pairs in semi-arid and arid regions (Georgescu et al 2011) (figure S3). In addition, urban-rural differences in background climate can exert an impact on vegetation growth (i.e. EVI values).

Thus, direct comparisons in EVI between urban and rural areas (e.g. urban EVI minus rural EVI) under drought conditions (PDSI < −2) cannot facilitate revealing differences in drought response of vegetation growth and productivity between urban and rural areas. Here, a linear mixed-effects model was used to help define drought response of vegetation growth and productivity since the model did not have a strong requirement for data distribution (i.e. independent variables may not be necessarily normally distributed, Schielzeth et al (2020)). First, we conducted a panel analysis to remove variation in EVI induced by weather variables including T_min, T_max, Prcp, and Srad using equation (1) for all urban and all rural areas individually (thus urban-rural difference in climate variables can be accounted for):

$$\begin{align}{y_{i,t}} &= \underbrace {{\alpha _1}t\! +\! {\mkern 1mu} {\alpha _2}T{\text{mi}}{{\text{n}}_{i,{} t}}\! +\! {\alpha _3}T{\text{ma}}{{\text{x}}_{i,{} t}} \!+\! {\alpha _4}{\text{Sra}}{{\text{d}}_{i,t}} \!+\! {\alpha _5}{\text{Prc}}{{\text{p}}_{i,t}}}_{{\text{Climatologically}}{} {\text{mean}}} \nonumber\\ &\quad+ \underbrace {(1|{\text{Cit}}{{\text{y}}_i})}_{{\text{EV}}{{\text{I}}_a}} + {\varepsilon _{i,t}}\end{align} \tag{ 1 }$$

where ${y_{i,t}}$ refers to EVI observations for city i in month t (June, July, or August) between 2001 and 2018, ${\alpha _1}$ captures the EVI trend over the study period among cities, ${\alpha _2},\,{\alpha _3},\,{\alpha _4},{\alpha _5}$ defines the sensitivity of EVI to T_min, T_max, Srad, and Prcp, respectively, $(1|{\text{Cit}}{{\text{y}}_i})$ accounts for the random effect without intercept among cities (each city as a group), and $\varepsilon$ is the error term. The component $(1|{\text{Cit}}{{\text{y}}_i})$ can account for variations that do not change with climate variables, e.g. variations in EVI attributed to spatial variations of soil quality and impervious surface fractions among different cities and rural counterparts, respectively. Although each urban-rural pair shares one PDSI value extracted from the monthly PDSI dataset, there may be still differences in drought conditions between urban and rural regions. The weather variables including T_min, T_max, and Precp in equation (1) can help account for such differences in drought conditions between urban and rural regions that may not be captured by the monthly PDSI dataset. Baseline models without weather variables or random effects were also tested; however, both Akaike information criterion and Bayesian information criterion values pointed to a better model as shown in equation (1).

In general, the grouped regression fitting through equation (1) provides a means to remove climate induced variation in EVI for each urban-rural pair and the model fitted values provided by ${\alpha _1} - {\alpha _5}$ terms are assumed as climatologically mean EVI values. Then, we define the difference between original EVI observations and climatologically mean EVI values (i.e. original minus mean) as EVI_a (figure 1 shows the flowchart deriving EVI_a). Thus, the drought resistance/response of vegetation growth and productivity (represented as Δx_a ) is defined using equation (2)

$$\begin{equation}\Delta {x_a} = {x_{a,d,\,l}} - {x_{a,nd,\,l}}.\end{equation} \tag{ 2 }$$

where x is EVI, ${x_{a,d}}$ refers to EVI_a under drought conditions (PDSI < −2), ${x_{a,nd}}$ stands for EVI_a under non-drought conditions (0 < PDSI < 2), and l represent either urban or rural landscapes.

With equation (2), any city-specific effects such as soil quality and impervious fractions on vegetation status can be largely removed, which provides more confidence in comparing vegetation growth and productivity among the selected 75 cities. Equation (2) was applied to urban and rural extents separately (75*2*3 = 450 times). Droughts are expected to lower the vegetation productivity, and thus $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ would be negative. A smaller reduction of EVI_a during major drought periods suggests a greater resistance of vegetation to the drought impact.

2.3. Urban-rural differences in environmental variables

To help explain the observed discrepancies in drought response of vegetation growth and productivity between urban and rural areas, variables associated with CO₂ and O₃ concentrations and mean temperature (T_m) were used. The mean temperatures for urban and rural extents individually were computed using T_min and T_max provided by the gridded Daymet product following equation (3):

$$\begin{equation}{T_{{\text{m}},{ }l}} = ({T_{{\text{max}},l}} + {T_{{\text{min}},l}})/2\end{equation} \tag{ 3 }$$

where l refers to either urban or rural extents. CO₂ concentrations within urban and rural extents were derived from the bias-corrected column-average dry air mole fraction of CO₂ (${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$) from Orbiting Carbon Observatory (Eldering et al 2017) (OCO-2, spatial resolution 1.3*2.25 km²). The ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ dataset was obtained from the reprocessed OCO-2 Lite files Version 10 r (available at https://disc.gsfc.nasa.gov/datasets) and has a retrieval accuracy of approximately 1 ppm. Although ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ is not a direct measurement of surface CO₂ concentration, the ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ dataset well characterizes the urban CO₂ dome (if any) at the city level (Kort et al 2012, Fu et al 2019) and urban-rural gradients in ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ have a linear relationship with surface CO₂ gradients (Wang et al 2019). Since the dataset is only available from 2014 and has large gaps in both spatial and temporal domains, the mean ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ values were computed for urban and rural extents individually using all observations available within the urban or rural extents (i.e. ${{\textrm X}_{{\text{C}}{{\text{O}}_2}}}$ values between 2015 and 2019). Thus, we did not further differentiate variation in urban-rural ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ gradients among JJA.

The ozone data were obtained using the hourly EPA's Air Quality Data (www.epa.gov/outdoor-air-quality-data) collected at outdoor monitors across the United States. Hourly O₃ data from each station within either an urban or rural polygon were aggregated to a monthly scale and then further averaged over all stations within that region (urban or rural). However, only around one-third of urban-rural pairs (24 of 75 urban-rural pairs) had valid monthly mean O₃ observations from 2001 to 2018, resulting in large data gaps for analysis. Thus, daily L2 total column O₃ data between 2004 and 2018 derived from the Ozone Monitoring Instrument (OMI) (Dobber et al 2006) onboard the Aura satellite (available at https://disc.gsfc.nasa.gov/datasets) were also used to facilitate surface-level O₃ estimations. The satellite based O₃ data were gridded at 0.25$^\circ$ by 0.25$^\circ$ and retrieved using the enhanced TOMS version-8 algorithm applied to the ultraviolet radiance data at 317.5 and 331.2 nm with a bias less than 3% (Balis et al 2007). We did not use the satellite based O₃ data directly for analysis in this study since the satellite O₃ dataset contained the total ozone column data. A regression analysis suggested that there was a statistically significant linear relationship between urban-rural differences in total column O₃ and surface O₃ gradients (equation (4)) (based on 24 urban-rural pairs, as shown in figure S5(C), with all three months of data considered). Thus, the regression equation was used to convert satellite based monthly mean O₃ differences between urban and rural areas to surface urban-rural differences in O₃. Behind this conversion, it was assumed that O₃ concentration was relatively stable over years for each month (either from satellite- or station-based observations). This assumption was evidenced by figures S5(A) and (B) showing that the ratio between standard deviation and mean of the O₃ concentration was relatively small (2%–10%) over years for each month within urban or rural extents.

The urban-rural differences in environmental variables (ΔE) including T_m, ${{\textrm X}_{{\text{C}}{{\text{O}}_2}}}$(CO₂), and O₃ were calculated using equation (4):

$$\begin{equation}{{\Delta E}} = {E_{\text{u}}} - {E_{\text{r}}}\end{equation} \tag{ 4 }$$

where E indicates the variable T_m, ${\text{C}}{{\text{O}}_2}$, or O₃, ${E_{\text{u}}}$ is the mean value in the urban extent, and ${E_{\text{r}}}$ refers to the mean value in the rural region. To determine the main factors in controlling differences in drought response of vegetation growth and productivity, the partial correlations of $\Delta {x_a}$ with latitude, longitude, urban size, mean monthly temperature, ΔT_m, $\Delta {\text{C}}{{\text{O}}_2}$, and ΔO₃ were computed. We binned urban-rural differences in T_m, ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}},$ and O₃ every 0.1 $^\circ {\text{C}}$, 0.1 ppm, and 0.1 ppbv to reduce stochastic error and data uncertainties. ΔT, ΔCO₂, ΔO₃ were also binned every 0.5 $^\circ {\text{C}}$, 0.5 ppm, and 0.5 ppbv; however, the new binning strategy would not change the significant linear slopes as shown in figures 4(A)–(C).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Drought resistance of vegetation growth and productivity in urban areas was stronger than drought resistance in rural areas based on analysis of monthly EVI from 2001 to 2018 (figure 2). Statistically significant differences (p-value < 0.05) were observed in $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ between urban and rural areas except for urban-rural pairs in early summer in arid climate (figure 2). For all selected urban-rural pairs, on average, $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}{ }$ in urban areas, was less negative than that in rural areas regardless of months within the growing season (JJA) (figure 2(B)). These negative $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ values are expected since $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ refers to the difference in ${\text{EV}}{{\text{I}}_{\text{a}}}$ between drought and non-drought conditions and mean ${\text{EV}}{{\text{I}}_{\text{a}}}$ under drought is smaller than that under non-drought conditions (figure S4(B)). A less negative value in $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ suggests a stronger ability to resist dampening effects of drought on vegetation growth. More specifically, mean $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ in urban areas was −0.0435, −0.0261, and −0.0174, and in rural areas, was −0.0779, −0.0535, and −0.0316, for JJA, respectively (figure 2(B)). Intuitively, it may be reasonable to assume that the less negative mean $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ value in urban areas stems from the fact that urban areas typically, except those in arid climate, exhibited a smaller EVI compared to rural regions in each month (JJA) (figure S3). The observed smaller EVI in urban areas in temperate and continental climate zones partly results from the coarse resolution of vegetated pixels (1 km) that suffer signal contaminations from underlying impervious surfaces (i.e. dampening effect of impervious surfaces on EVI). This mixed signal effects (or spectral mixture of vegetation and impervious surfaces within a pixel) can further propagate to ${\text{EV}}{{\text{I}}_{\text{a}}}$ and $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ calculations, causing smaller standard deviations in ${\text{EV}}{{\text{I}}_{\text{a}}}$ (figure S4) and $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ in urban regions compared to rural regions (figure 2). However, such an intuitive explanation should not be a concern since the $\Delta {\text{EV}}{{\text{I}}_{\text{a}}}$ was computed as the difference in ${\text{EV}}{{\text{I}}_{\text{a}}}$ between drought and non-drought conditions for urban and rural extents individually (see the section 2 for further details), i.e. we did not compare ${\text{EV}}{{\text{I}}_{\text{a}}}$ under drought conditions in urban areas with ${\text{EV}}{{\text{I}}_{\text{a}}}$ under drought conditions in rural areas directly.

As regional climate largely controls the types of ecosystems, we further grouped urban-rural pairs into three main climate zones based on Köppen-Geiger scheme (Beck et al 2018), including arid, temperate, and continental climates. The drought resistance of vegetation growth within each climate zone is shown in figures 2(C)–(E). The separate analysis among climate zones and months shows consistent findings that urban vegetation exhibits stronger ability to resist drought impacts than rural vegetation though the urban-rural differences in vegetation drought resistance varies (figures 2(C)–(E)). In the arid zone, on average, ${\text{dEV}}{{\text{I}}_{\text{a}}}$ (urban-rural differences in drought resistance of vegetation growth) was 0.0047 and 0.0135 for July and August (p-value < 0.01), respectively while there was no difference in ${\text{dEV}}{{\text{I}}_{\text{a}}}$ for June (p-value > 0.01). Compared to the arid zone, the temperate and continental climate zones in general showed a much higher value for ${\text{dEV}}{{\text{I}}_{\text{a}}}$ in June and July. For example, the mean ${\text{dEV}}{{\text{I}}_{\text{a}}}$ in July for urban-rural pairs under the continental climate was 0.0133 and under the temperate climate was 0.0177, much higher than 0.0047 for urban-rural pairs from the arid zone. For mean ${\text{dEV}}{{\text{I}}_{\text{a}}}$ in August, urban-rural pairs from the temperate zone exhibited the largest value (0.0148), followed by those from the arid (0.0135) and continental (0.017) climate zones.

We next explored the potential effects of several drivers underlying urban-rural discrepancies in drought resistance of vegetation growth, including urban-rural differences in monthly mean temperature (ΔT), CO₂ (ΔCO₂), and O₃ (ΔO₃), city size (Size), the latitude and longitude of a city (Lat and Lon, in case spatial location would be a factor), and background monthly mean temperature (T_m). To find the main factors controlling ${\text{dEV}}{{\text{I}}_{\text{a}}}$, partial correlation coefficients (partial CC) between ${\text{dEV}}{{\text{I}}_{\text{a}}}$ and potential drivers were computed for all selected urban-rural pairs from all three climate zones (figure 3(A)) as well as separately for arid (figure 3(B)), temperate (figure 3(C)), and continental climates (figure 3(D)). The main factors in controlling ${\text{dEV}}{{\text{I}}_{\text{a}}}$ across the conterminous U.S. were ΔT (partial CC of 0.36), ΔCO₂ (0.31), and ΔO₃ (−0.12) as seen in figure 3(A). Significant linear relations were also observed for the associations of ΔT, ΔCO₂, ΔO₃ and ${\text{dEV}}{{\text{I}}_{\text{a}}}$ from binned observations (figure 4) that were adopted to reduce stochastic error and data uncertainties in the analysis. The main controlling factors for ${\text{dEV}}{{\text{I}}_{\text{a}}}$ in each climate zone were different among three climate backgrounds. For example, the most influential variable was ΔO₃ in the arid zone (partial CC of 0.37, figure 3(B)), while in the temperate zone, the main drivers influencing ${\text{dEV}}{{\text{I}}_{\text{a}}}$ were ΔCO₂ (partial CC of 0.22), ΔO₃ (partial CC of −0.18), and T_m (partial CC of −0.15) (figure 3(C)). In the continental zone (figure 3(D)), ΔCO₂, followed by ΔO₃ was the main factor in affecting ${\text{dEV}}{{\text{I}}_{\text{a}}}$. Overall, ΔT was identified as a main driver for ${\text{dEV}}{{\text{I}}_{\text{a}}}$ over the CONUS but not in each climate zone. Further explanations for this finding were presented in the discussion section.

We found urban ecosystems were more resistant to drought compared to their rural counterparts regardless of climate zones across CONUS. However, the main variables that contribute to such urban-rural differences in drought resistance of vegetation growth varies among climate zones. In general, surface temperature gradient (ΔT) was the most influential factor (p-value < 0.01, figure 2(A)) at the continental scale that consists of various ecosystems. Higher temperature in urban areas in relative to rural areas, observed for most of the selected urban-rural pairs (figure S6), likely leads to an earlier start but later end of the growing season and thus a longer growing season (Li et al 2017, Wang et al 2019), resulting in enhanced vegetation growth prevalent in urban areas (Zhao et al 2016) (while still a lower EVI value compared to rural areas due to the dampening effect from impervious surfaces). Thus, when the very same drought occurred in both urban and rural areas, the enhanced vegetation growth attributed to higher temperatures in urban areas exhibited a better ability/chance to endure/survive drought stress. In addition, plants growing in warmer/drier urban environments may already have had (screened before planting) or developed traits that allow them to better cope with drought stress (Anderegg and HilleRisLambers 2016). However, ΔT was not the controlling factor driving the urban-rural contrast of vegetation drought resistance within individual climate zones as shown in figure 3 partly due to the much less variation of ΔT in each climate zone (figure S7). Specifically, as shown in figure S7, vegetation growth in arid or semi-arid cities in general exhibits a much stronger drought resistance compared to that in the other two climate zones despite a wide range of ΔT from −2.0 °C to 3.5 $^\circ {\text{C}}$. Even though urban-rural difference in drought resistance varies from −0.10 to 0.11 within temperate or continental climate regions, in a similar magnitude of variation in arid climate, there is a much smaller variation of ΔT. Thus, the urban-rural temperature gradient is not strong enough to cause such considerable differences in vegetation growth. The stronger drought resistance of vegetation (mainly shrub and grass) growth in cities from the arid zone was mainly related to vegetation growth in rural areas exposed to higher ozone concentration (figure S8). Higher ozone concentration in rural areas than in cities (Gregg et al 2003), could cause more reduction in photosynthesis and vegetation growth in rural areas (Ainsworth et al 2012), which accounted for the negative relationship between surface O₃ gradients and urban-rural differences in drought resistance of vegetation growth (figures 3 and 4(C)). In addition, urban-rural difference in drought resistance of vegetation growth in the temperate zone is also sensitive to background monthly mean temperature T_m. Since the selected urban-rural pairs from the temperate climate zone exhibits the highest background monthly mean temperature (figure S9), this result emphasizes possible temperature stress on vegetation growth, particularly in urban areas, given relatively stable temperature gradients for urban-rural pairs from the temperate zone (figure S7).

Surface CO₂ gradient is identified as another important factor, in addition to surface temperature and O₃ gradients, responsible for the observed urban-rural discrepancies in drought resistance of vegetation growth. The surface CO₂ gradient is represented by the satellite-based ${{\textrm X}_{{\text{C}}{{\text{O}}_{\text{2}}}}}$ gradient since there is a positive, linear relationship between the two (with a scale factor of ∼25) (Wang et al 2019). This satellite-based CO₂ gradient dataset (figure S10) shows that a greater difference in CO₂ concentration between urban and rural areas could result in a much stronger urban-rural difference in drought resistance of vegetation growth, particularly under temperate and continental climates (figures 3(A), (C), (D) and 4(B)). This conclusion is consistent with previous studies that emphasize the beneficial effect of atmospheric CO₂ enhancement on vegetation growth under drought conditions by stimulating photosynthesis (Ainsworth and Rogers 2007) or by increasing the water use efficiency of plants (Keenan et al 2013). However, the beneficial effect of CO₂ on urban-rural contrast in drought resistance is observed in temperate and continental zones rather than in arid zones even though urban-rural pairs from the arid zone typically showed the largest CO₂ gradient (figure S11). The insensitivity of drought impact on vegetation growth to CO₂ gradients in the arid zone, relative to other two climate zones, may result from other more limiting factors such as nutrient deficiency (Wang et al 2020) or stomatal closure in response to rising CO₂ at the cost of enhanced growth (Frank et al 2015) in urban regions from the arid zone. Drought in the water-limited region is often more severe (as indicated in figure S12), and it is also possible that such drought condition completely overwhelms the beneficial effects of elevated CO₂ on vegetation growth.

Despite the identified main variables, factors such as landscape configuration/composition, atmospheric deposition (e.g. nitrogen and phosphorus), and management practices may also affect the observed discrepancies in drought resistance of vegetation growth between urban and rural areas. For example, urban landscape configuration and composition has been related to surface temperature gradients between urban and rural areas (Connors et al 2013, Estoque et al 2017), thus either strengthening or accentuating the effect of temperature gradients on difference in drought resistance between urban and rural vegetation growth. As urban areas typically have a higher atmospheric nitrogen and phosphorus deposition (Decina et al 2018), this may lessen nutrient restrictions on vegetation growth, thus contributing positively to the observed differences in drought resistance of vegetation growth between urban and rural areas. Human practices such as urban irrigation, however, may be a factor weakening the conclusion made towards the main factors driving the urban-rural differences in drought resistance of vegetation growth. Thus, we repeated the analysis by excluding some urban-rural pairs from the arid zone where irrigation typically was most pronounced and led to 'the oasis effect' (Georgescu et al 2011), i.e. higher EVI in cities than in rural areas (figure S3). Further analysis of urban-rural pairs from the three climate zones suggests that the main factors controlling the urban-rural differences in drought resistance of vegetation growth are still gradients of temperature and CO₂ and O₃ concentrations but with a smaller partial CC (figures S13 and 3(A)). Thus, although urban irrigation contributed positively to the observed urban-rural contrasts in drought resistance of vegetation growth, such a factor alone could not explain all the observed variances associated with drought resistance. Even under strong irrigation for cities in the arid zone, surface O₃ gradient was still identified as the main driver for observed discrepancies in drought resistance of vegetation growth between urban and rural areas (figure 3(B)). As such, there is strong evidence in attributing the observed discrepancies in drought resistance of vegetation growth between urban and rural areas to factors including temperature, CO₂, and O₃ gradients. Urban-rural gradients in wind speed and humidity were also accounted for in the analysis (figure S15). The results suggested the controlling factor for the urban-rural differences in drought resistance would still be temperature, CO₂, and O₃ gradients (figure S15).

Benefiting from long-term satellite and station-based datasets, our study provides a data driven approach for revealing the impacts of competing and interactive environmental factors on response of vegetation growth to drought at a large scale. This approach advocates the urban plant physiology concept by utilizing urban-rural contrasts in environmental conditions, such as altered temperature gradients and O₃ concentration. The substantial differences in environmental conditions along urban-rural gradients can provide a 'natural laboratory' that is more widely accessible, compared to manipulative environments, to the scientific community to understand ecosystem traits under a changing climate. With projected increase in temperature and atmospheric CO₂ in the future, we can also select only urban-rural pairs with positive values in ΔT and ΔCO₂ for analysis. In this case, it is found that the mean urban-rural difference in drought resistance of vegetation growth becomes positive, i.e. 0.0389 for June, 0.0294 for July, and 0.0201 for August (figure S14). Thus, the novel aspect of this study is to provide a more general assessment of vegetation responses to environmental factors across different regions, and ultimately, to reduce uncertainties in quantifying terrestrial carbon dynamics.

Fu and Bernacchi would like to acknowledge the funding support from Global Change and Photosynthesis Research Unit of the USDA Agricultural Research Service. Hu is partly funded by NASA Projects (80NSSC20K1263 and 80NSSC21K0430). Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the author(s) and do not necessarily reflect the views of the USDA. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture. USDA is an equal opportunity provider and employer. We would like to thank Jesse McGrath for writing R scripts in accessing US EPA Air Quality Data. We are also grateful to editors and three anonymous reviewers for their thought comments and suggestions which helped improve this manuscript.

All data that support the findings of this study are included within the article (and any supplementary files).

The authors declare no conflicts of interest relevant to this study.