This is the accepted version of the article: Wang, Songhan; Ju, Weimin; Peñuelas, Josep; [et al.]. «Urban−rural gradients reveal joint control of elevated CO2 and temperature on extended photosynthetic seasons». Nature ecology and evolution, Vol. 3 (March 2019), p. 1076–1085. DOI 10.1038/s41559-019-0931-1 This version is avaible at https://ddd.uab.cat/record/218190 under the terms of the license Supplementary Information for Urban−rural gradients reveal joint control of elevated CO2 and temperature on extended photosynthetic seasons Songhan Wang, Weimin Ju, Josep Peñuelas, Alessandro Cescatti, Yuyu Zhou, Yongshuo Fu, Alfredo Huete, Min Liu, Yongguang Zhang Yongguang Zhang Email: yongguang_zhang@nju.edu.cn This PDF file includes: Supplementary Figure 1 to 18 Supplementary Table 1 to 3 References for SI reference citations 1 Supplementary Figure 1 | Locations of the study areas and background land-cover types. Red dots indicate the 880 urban clusters at mid–high latitudes in the Northern Hemisphere (≥30°). The background colours represent the land-cover types extracted from the Climate Change Initiative product for 2015. EBF: evergreen broadleaf forest; DBF: deciduous broadleaf forest; ENF: evergreen needleleaf forest; DNF: deciduous needleleaf forest. 2 Supplementary Figure 2 | Seasonal variation of SIF and EVI in urban and rural areas. a, Normalized monthly EVI and SIF seasonal cycles for urban and rural areas from 2015 to 2017 at mid–high latitudes in the Northern Hemisphere. We chose 30-km rural buffers as the representative area for rural conditions. Error bars indicate monthly standard deviations. p indicates the significance of the difference between urban and rural SIF or EVI (two-sample Kolmogorov–Smirnov test). b, Monthly normalized urban–rural ∆SIF and ∆EVI (∆SIF = SIFurban − SIFrural and ∆EVI = EVIurban − EVIrural). The values in urban areas are higher than those in rural areas during spring and autumn for both SIF and EVI. The urban–rural differences of SIF are significantly higher than those of EVI. 3 4 Supplementary Figure 3 | SIF- and EVI-based phenological indicators (SOS, POS, EOS and GSL) for all urban clusters and their corresponding rural buffers (n = 880). In each plot, the x-axis shows the timings of the relevant phenology indicator in urban areas, and the y-axis shows the timings in the surrounding rural areas (10-40 km buffers, at distances of 10, 20, 30 and 40 km from the relevant corresponding urban area). The columns represent the phenological indicators (SOS: start of growing season; POS: peak of growing season; EOS: end of growing season; GSL: length of growing season). The black dotted lines represent the 1:1 lines, and the solid red lines represent the regression lines between urban and rural phenological indicators. The slopes and 95% confidence intervals of each regression line are shown. Slopes greater than 1 indicate that the phenological stage is occurring at a later date in rural areas than urban areas. ∆ represents the mean difference between urban areas and rural buffers. In all cases, p < 0.001 (p indicates the significance of the linear regression). 5 Supplementary Figure 4 | Mean differences in phenological metrics for SIF and EVI between urban areas and corresponding rural buffer zones of increasing distances from the city. The operator ∆ represents the differences in phenological indicators between urban areas and rural buffers (e.g., ∆SOS = SOSurban − SOSrural). Error bars represent standard errors of the mean (n = 880); p represents the significance between the mean value and zero; *: p < 0.05. SOS: start of growing season; POS: peak of growing season; EOS: end of growing season; GSL: length of growing season. 6 Supplementary Figure 5 | Mean urban–rural phenological differences based on SIF and EVI, which have the same 5-km spatial resolution. We firstly resampled the OCO-2 SIF and MODIS EVI to a spatial resolution of 5 km, and next calculated the phenology indicators in urban clusters and corresponding rural buffers following the same procedure in the main text. SOS, start of growing season; POS, peak of growing season; EOS, end of growing season; GSL, length of growing season. The operator ∆ represents the differences in phenological indicators between urban areas and rural buffers (e.g., ∆SOS = SOSurban − SOSrural). The error bars represent standard errors of the mean (n = 880). 7 Supplementary Figure 6 | Mean urban–rural phenological differences (∆SOS, ∆POS, ∆EOS and ∆GSL) across different climate zones. The x-axis represents different climate zones (i.e., arid, temperate and boreal) from the Köppen–Geiger classification zones. The y-axis represents the mean urban–rural phenological differences. Error bars represent standard errors of the mean. SOS: start of growing season; POS: peak of growing season; EOS: end of growing season; GSL: length of growing season. 8 Supplementary Figure 7 | Land surface temperature (LST) and CO2 concentration changes from 2015 for four projected Representative Concentration Pathway (RCP) scenarios. a, Mean LST changes from 2015 at mid- to high-latitudes in the Northern Hemisphere (≥30°) simulated by four Earth system models (CSIRO-Mk3.6.0, GFDL-CM3, GISS-E2-H and NorESM1-M) from the Coupled Model Intercomparison Project Phase 5 (CMIP5). b, Projected CO2 concentration changes from 2015. 9 Supplementary Figure 8 | Urban–rural phenological differences (∆GSL) across different environmental factors: (a) precipitation; (b) urban–rural NO2 difference (∆NO2); (c) urban–rural O3 difference (∆O3). Based on precipitation, ∆NO2 and ∆O3, the 880 urban clusters are classified into four categories using different quartiles. The mean values and standard errors of the mean are shown with bars and error-bars. The dashed lines represent the linear regression of the mean values. ∆NO2 = NO2 urban – NO2 rural; ∆O3 = O3 urban – O3 rural. GSL: length of growing season. 10 Supplementary Figure 9 | Results after excluding the urban clusters in arid climate zones. a, Partial correlations of SIF urban–rural phenological differences (∆SOS, ∆POS, ∆EOS and ∆GSL) with ∆Tday and ∆CO2, based on all the 880 urban clusters (‘Original’) and the remaining urban clusters after excluding those in the arid climate zones (‘Filtered’). ∆Tday: daytime urban–rural land surface temperature difference; ∆CO2: urban–rural CO2 difference. SOS: start of growing season; POS: peak of growing season; EOS: end of growing season; GSL: length of growing season. *: p < 0.05; **: p < 0.01. b, Projected shifts of GSL in per decade (days·decade-1). The error bars represent standard errors of the mean. 11 Supplementary Figure 10 | Partial correlations of SIF urban–rural30km and rural10km–rural30km phenological differences with ∆Tday and ∆CO2. Urban, rural10km and rural30km represent the phenology indicators in the urban clusters, 10-km rural buffers and 30-km rural buffers, respectively. ∆Tday: daytime urban–rural land surface temperature difference; ∆CO2: urban–rural CO2 difference. SOS: start of growing season; POS: peak of growing season; EOS: end of growing season; GSL: length of growing season. *: p < 0.05; **: p < 0.01 (n = 880). 12 Supplementary Figure 11 | An example of the urban cluster and its corresponding rural buffers using Moscow. a, The urban cluster and rural buffers (10-, 20-, 30- and 40- km) and the background land cover classifications. The urban cluster was derived from the satellite nighttime light data using a cluster-based method, which defined the boundary between urban and rural areas, regardless of the internal spatial heterogeneity. The rural buffers were defined by extending 10, 20, 30 and 40 km outward from the urban cluster perimeter. We excluded the pixels that were water bodies or crops in this analysis. b, The pixels excluded using the elevation mask. We also excluded the pixels in the rural areas which had elevations greater than ±50 m of the average elevation of urban pixels. 13 Supplementary Figure 12 | Relative retrieval uncertainties of single SIF measurements in urban clusters and corresponding rural buffers. The relative retrieval uncertainties were calculated as the ratios between single-measurement precision and the absolute SIF values, which were derived from the original OCO-2 SIF products. The mean values and standard errors (SEs) of the mean are shown with bars and error-bars. 14 Supplementary Figure 13 | R2 values between the observation time series and double-sigmoidal fits for SIF (a) and EVI (b). Most of the R2 values for SIF are greater than 0.8. Only a few urban clusters have R2 values between 0.75 and 0.8. Most of the R2 values for EVI are greater than 0.9, which are higher than those for SIF. 15 Supplementary Figure 14 | Robustness of the phenology extraction method. a, An example of the extraction of the phenological indicators. The gray points represent the SOS, POS and EOS respectively. SOS, start of growing season; POS, peak of growing season; EOS, end of growing season; GSL, length of growing season. b, Extraction of the phenological indicators by changing the moving window, i.e., we assigned a weight of two if a central point was within ±50% of the median for a moving window of five points. c, Extraction of the phenological indicators without using the weighting scheme. d, Similar with Fig. 1 in the main text, except the results were calculated by setting the moving window to five points. The operator ∆ represents the differences in phenological indicators between urban areas and rural buffers (e.g., ∆SOS = SOSurban − SOSrural). The error bars represent standard errors of the mean (n = 880). e, Similar with Fig. 1 in the main text, except the results were calculated without using the weighting scheme. 16 Supplementary Figure 15 | Spatial distribution of near-surface CO2 concentrations in Shanghai from 172 sampling points. 17 Supplementary Figure 16 | Factors for converting urban–rural XCO2 gradients into surface CO2 gradients. a, Atmospheric XCO2 observations across Shanghai as a function of the distance to the city centre. b, Surface CO2 measurements in Shanghai as a function of the distance to the city centre. The factor for converting XCO2 gradients into surface CO2 gradients (K) is SlopeCO2/SlopeXCO2 = 21.2. c, Scatter plots between urban–rural observed XCO2 gradients and surface CO2 gradients from previous studies, with K = 28.5. 18 Supplementary Figure 17 | Partial correlations of urban–rural EOS gradients with CO2 gradients (a) and predicted GSL shifts based on RCP 4.5 scenario (b) using different K values. K represents the factor for converting XCO2 gradients into surface CO2 gradients. The K factor derived from the first approach was 21.2 ± 8.3; while that from the second approach was 28.5 ±9.5 (mean ±standard error). Therefore we used a series of K factors ranging from 19 to 30, which was the intersection of the K factors estimated from these two approaches. ∆CO2: urban–rural CO2 difference. EOS: end of growing season; GSL: length of growing season. *: p < 0.05 (n = 880). The error bars in panel b represent the standard errors (SEs) of the mean. 19 Supplementary Figure 18 | Comparisons of the future GSL (length of growing season) shifts based on urban–rural land surface temperature (LST) gradients and air surface temperature (AST) gradients. a, Mean urban–rural LST difference and AST difference. b, sensitivities of GSL to LST and AST. c, Projected shifts of GSL in per decade (days·decade-1) based on urban–rural LST differences and AST differences, respectively. The error bars represent standard errors of the mean. 20 Supplementary Table 1 | Descriptions of the data sets. Variable Origin Temporal Spatial Date Reference resolution resolution Urban clusters – – 1 km 2013 Zhou, et al. 1 SIF (solar-induced OCO-2 (version 8r) Approximately 1.3×2.25 km2 2015–2017 Frankenberg, et al. 2 chlorophyll 16 days fluorescence) EVI (enhanced MODIS (MYD13A2) 16 days 1 km 2015–2017 Huete, et al. 3 vegetation index) LST (land surface MODIS (MYD11A2) 16 days 1 km 2015–2017 Wan and Dozier 4 temperature) XCO2 OCO-2 (version 8r) approximately 1.3×2.25 km2 2015–2017 Crisp, et al. 5 (column-averaged CO2 16 days mixing ratios) Air temperature GHCN-CAMS monthly 0.5° 2015–2017 Fan and Van den Dool 6 analyzed data Precipitation GPM-IMERG data monthly 0.1° 2015–2017 Huffman, et al. 7 Land cover ESA CCI projects – 300 m 2015 ESA 8 Altitude GTOPO30 DEM – 1 km 1996 Gesch, et al. 9 AST (air surface TerraClimate monthly 4 km 2015–2017 Abatzoglou, et al. 10 temperature) NO2 Ozone Monitoring monthly 0.125° 2015–2017 Boersma, et al. 11 Instrument (OMI) O3 Ozone Monitoring monthly 0.25° 2015–2017 Veefkind, et al. 12 Instrument (OMI) 21 Supplementary Table 2 | Observations of urban–rural surface CO2 gradients from previous studies. The urban–rural gradient for Phoenix ranges from 19 to 185 ppm, which is highly variable. Thus, we excluded these observations when calculating the conversion factor. The urban–rural gradients from these studies averaged approximately 30 ppm (excluding the Phoenix data). No. City, country Year Observation method Urban CO2 Rural CO2 Urban–rural Relative CO2 Reference concentration (ppm) concentration (ppm) CO2 gradient (ppm) increase (%) 1 Phoenix, USA 1998 Mobile measurements 555 370 185 50.0 Idso, et al. 13 2 Phoenix, USA 2000 Mobile measurements 619.3(max) 424.3–490.6 – 41.5 Idso, et al. 14 3 Phoenix, USA 2002 Site observations 396 377 19 5.0 Day, et al. 15 4 Essen, Germany 2002–2003 Mobile measurements 404 371 33 8.9 Henninger and Kuttler 16 5 Essen, Germany 2010 Site observations 417 381 36 9.4 Büns and Kuttler 17 6 Rome, Italy 1995–2004 Different zones 477 405 72 17.8 Gratani and Varone 18 7 Baltimore, USA 2007 Site observations 488 422 66 15.6 George, et al. 19 8 Valladolid, Spain 2006–2007 Site observations 397.8 393.3 4.5 1.1 Garcí a, et al. 20 9 Portland, USA 2009 Site observations 408 403 5 1.2 Rice and Bostrom 21 10 Boston, USA 2011 Site observations 414 393 21 5.3 Briber, et al. 22 11 Basel, Switzerland 2005–2015 Site observations – – 10.5–12 – Schmutz, et al. 23 12 Berkeley, USA 2013 Site observations 425 407 18 4.4 Shusterman, et al. 24 13 Nanjing, China 2008 Mobile measurements 435 387 48 12.4 Zhao, et al. 25 14 Xiamen, China 2011 Mobile measurements 477.3 434.1 43.2 10.0 Li, et al. 26 15 Shanghai, China 2014 Mobile measurements 469.2 414.1 55.1 13.3 Zhu, et al. 27 22 Supplementary Table 3 | Urban–rural photosynthetic phenological differences (SIF) and explanatory variables based on forward stepwise regression models. SOS: start of season; POS: peak of season; EOS: end of season; GSL: length of season; ∆LSTday: daytime urban–rural land surface temperature difference; ∆CO2: urban–rural CO2 difference; SE: standard error. SIF Model Explanatory variables Coefficient (±SE) R2 p ∆SOS 1 ∆LSTday (°C) –1.54 ±0.60 0.11 0.008 2 ∆LSTday (°C) –1.55 ±0.60 0.15 0.004 ∆CO2 (ppm) –0.15 ±0.07 ∆POS 1 ∆LSTday (°C) –2.10 ±0.79 0.12 0.008 2 ∆LSTday (°C) –2.11 ±0.78 0.15 0.004 ∆CO2 (ppm) –0.18 ±0.09 ∆EOS 1 ∆CO2 (ppm) 0.17 ±0.07 0.10 0.023 ∆GSL 1 ∆LSTday (°C) 3.33 ±1.11 0.13 0.003 2 ∆LSTday (°C) 3.31 ±1.10 0.19 0.000 ∆CO2 (ppm) 0.32 ±0.11 23 References 1 Zhou, Y. et al. A global map of urban extent from nightlights. Environmental Research Letters 10, 054011 (2015). 2 Frankenberg, C. et al. Prospects for chlorophyll fluorescence remote sensing from the Orbiting Carbon Observatory-2. Remote Sensing of Environment 147, 1-12 (2014). 3 Huete, A. et al. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote sensing of environment 83, 195-213 (2002). 4 Wan, Z. & Dozier, J. A generalized split-window algorithm for retrieving land-surface temperature from space. IEEE Transactions on geoscience and remote sensing 34, 892-905 (1996). 5 Crisp, D. et al. The on-orbit performance of the Orbiting Carbon Observatory-2 (OCO-2) instrument and its radiometrically calibrated products. Atmospheric Measurement Techniques 10, 59 (2017). 6 Fan, Y. & Van den Dool, H. A global monthly land surface air temperature analysis for 1948–present. Journal of Geophysical Research: Atmospheres 113 (2008). 7 Huffman, G. J. et al. NASA global precipitation measurement (GPM) integrated multi-satellite retrievals for GPM (IMERG). Algorithm theoretical basis document, version 4, 30 (2015). 8 ESA. Land Cover CCI Product User Guide Version 2.0. (2017). 9 Gesch, D. B., Verdin, K. L. & Greenlee, S. K. New land surface digital elevation model covers the Earth. EOS, Transactions American Geophysical Union 80, 69-70 (1999). 10 Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Scientific data 5, 170191 (2018). 11 Boersma, K. et al. Near-real time retrieval of tropospheric NO2 from OMI. Atmospheric Chemistry and Physics 7, 2103-2118 (2007). 12 Veefkind, J. P., de Haan, J. F., Brinksma, E. J., Kroon, M. & Levelt, P. F. Total ozone from the Ozone Monitoring Instrument (OMI) using the DOAS technique. IEEE Transactions on Geoscience and Remote Sensing 44, 1239-1244 (2006). 13 Idso, C. D., Idso, S. B. & Balling Jr, R. C. The urban CO 2 dome of Phoenix, Arizona. Physical Geography 19, 95-108 (1998). 14 Idso, C. D., Idso, S. B. & Balling Jr, R. C. An intensive two-week study of an urban CO2 dome in Phoenix, Arizona, USA. Atmospheric Environment 35, 995-1000 (2001). 15 Day, T. A., Gober, P., Xiong, F. S. & Wentz, E. A. Temporal patterns in near-surface CO2 concentrations over contrasting vegetation types in the Phoenix metropolitan area. Agricultural and Forest Meteorology 110, 229-245 (2002). 16 Henninger, S. & Kuttler, W. Near surface carbon dioxide within the urban area of Essen, Germany. Physics and Chemistry of the Earth, Parts A/B/C 35, 76-84 (2010). 17 Büns, C. & Kuttler, W. Path-integrated measurements of carbon dioxide in the urban canopy layer. Atmospheric environment 46, 237-247 (2012). 18 Gratani, L. & Varone, L. Daily and seasonal variation of CO2 in the city of Rome in relationship with the traffic volume. Atmospheric Environment 39, 2619-2624 (2005). 19 George, K., Ziska, L. H., Bunce, J. A. & Quebedeaux, B. Elevated atmospheric CO 2 concentration 24 and temperature across an urban–rural transect. Atmospheric Environment 41, 7654-7665 (2007). 20 García, M. Á., Sánchez, M. L. & Pérez, I. A. Differences between carbon dioxide levels over suburban and rural sites in Northern Spain. Environmental Science and Pollution Research 19, 432-439 (2012). 21 Rice, A. & Bostrom, G. Measurements of carbon dioxide in an Oregon metropolitan region. Atmospheric Environment 45, 1138-1144 (2011). 22 Briber, B. M., Hutyra, L. R., Dunn, A. L., Raciti, S. M. & Munger, J. W. Variations in atmospheric CO2 mixing ratios across a Boston, MA urban to rural gradient. Land 2, 304-327 (2013). 23 Schmutz, M., Vogt, R. & Parlow, E. in 9th International Conference on Urban Climate 101-105 (Toulouse, France, 2015). 24 Shusterman, A. A. et al. The BErkeley Atmospheric CO2 Observation Network: initial evaluation. Atmospheric Chemistry and Physics 16, 13449-13463 (2016). 25 Zhao, H., Ouyang, y. & Qi, j. Near-surface atmospheric properties along an urban center-suburban-urban forest gradient in summer in Nanjing City. ACTA ECOLOGICA SINICA 29, 6654-6663 (2009). 26 Li, Y., Xing, Z. & Mu, C. Spatial and Temporal Variations of Near Surface Atmospheric CO 2 with Mobile Measurements in Fall and Spring in Xiamen, China. ENVIRONMENT SCIENCE 35, 1671-1679 (2014). 27 Zhu, X., Pan, C. & Liu, M. Spatial characristics of near surface CO 2 concentrations and analysis of its influencing factors in spring in Shanghai city. Resources and Environment in the Yangtze Basin 24 (2015). 25