Projecting the CO 2 and Climatic Change Effects on the Net Primary Productivity of the Urban Ecosystems in Phoenix , AZ in the 21 st Century under Multiple RCP ( Representative Concentration Pathway ) Scenarios

Urban vegetation provides ecological services that promote both the ecosystem integrity and human well-being of urban areas, and thus is critical to urban sustainability. As a key indicator of ecological health, net primary productivity (NPP) provides valuable information about the performance of urban ecosystem in response to the changes in urban climate and atmosphere in the 21st century. In this study, a process-based urban ecosystem model, HPM-UEM (Hierarchical Patch Mosaic-Urban Ecosystem Model), was used to investigate spatiotemporal dynamics of urban ecosystem NPP in the Phoenix city, AZ under three representative concentration pathway (RCP2.6, RCP4.5 and RCP8.5) during the 21st century. The results indicated that, by the end of the 21st century, the urban ecosystem’s NPP would increase by 14% (in RCP2.6), 51% (in RCP4.5) and 99% (in RCP8.5) relative to that in the late 2000s, respectively. Factorial analysis indicated that CO2 fertilization effect would be the major driver of NPP change, accounting for 56–61% of the NPP increase under the scenarios. Under the RCP2.6 scenario, the strongest NPP increase would be found in the agricultural lands located in the west and southeast of the city. Under the RCP4.5 and RCP8.5 scenarios, the strongest NPP increase would be found in the mesic residential areas that mainly located to the eastern, southern, and southwestern of the Phoenix Mountains Preserve. Although higher ecosystem NPP in the future implies improved ecosystem services that may help to alleviate the heat stress (by providing more shading) and air pollution in the city, this will be at the cost of higher irrigation water usage, probably leading to water shortage in the natural ecosystems in this arid region. Furthermore, this study indicated the rich (such as in mesic residential area) would enjoy more benefits from the improved urban ecosystem services than the poor (such as in xeric residential area).


Introduction
Since the industrial revolution, the earth has entered an era in which human impacts on the fluxes of radiative energy through the earth system have demonstrable effects on climate (Intergovernmental Panel on Climate Change (IPCC) 2007).Net primary productivity (NPP), or the rate of assimilation of carbon dioxide (CO 2 ) through photosynthesis, is the fundamental link between the atmosphere and the biosphere [1].NPP contributes to human survival because it is not only the basis for food, fiber and wood production, but also affects composition of the atmosphere, the availability of fresh water, biodiversity, and the adjusting mechanisms of energy supply and distribution [2].NPP is a sensitive indicator of climatic and other forms of environmental change [3], and is increasingly pertinent to land-use management [2,4].Vackar et al. found that the human appropriation of NPP (HANPP, which is an important indicator of environmental sustainability) amounted to 56% of the annual potential natural NPP in the territory of the Czech Republic [5].In particular, urban ecosystem NPP is a key indicator of vegetation health in urban ecosystems, which provide critical services to urban residents, including the "regulating services" (e.g., purification of air and water, and regulation of climate and floods), "cultural services" (e.g., recreational, spiritual, religious and other nonmaterial benefits), and "supporting services" (e.g., nutrient cycling) [6].For example, Lazzarini et al. found the observed diurnal cool island effect in hot desert cities could be largely explained by relative vegetation abundance [7].Since human well-being in urban areas depends on the ecosystem services provided by urban vegetation, urban sustainability depends on the sustainability of urban ecosystems [8].
The ecosystem health of urban vegetation and the ecosystem services provided by them may be vulnerable because of the intensified environmental stresses (e.g., heat and drought) due to future climate change [9].Multiple lines of evidence indicated that urban ecosystem NPP is sensitive to environmental changes at various spatiotemporal scales.Previous studies have investigated inter-annual variability of NPP in response to climate change [10][11][12][13], the contributions of urbanization and climate change to NPP variations [14][15][16][17], and the effect of urban land transformation on NPP [18].Most research has been focused on the effect of land use/land cover change and historical climate change; however, the dynamic of NPP in response to future climate/CO 2 change effects are still unclear.To predict and manage urban ecosystem sustainability in the 21st century, it is important to estimate the changes of urban ecosystem NPP in response to various future trajectories of climate and CO 2 .Ecosystem process-based modeling is a powerful tool to quantitatively estimate and predict dynamics in terrestrial biogeochemical cycles and isolate/analyze associated control mechanisms in the context of global change [19,20].
In this study, an urban ecosystem model, the Hierarchical Patch Mosaic-Urban Ecosystem Model (HPM-UEM), was used to predict the responses of the NPP of the Phoenix city, AZ, USA to climate and CO 2 changes in the 21st century.The model simulations were conducted in 500-m spatial resolution, and the NPP changes under three representative concentration pathways-RCP2.6,RCP4.5 and RCP8.5-were compared.Numeric experiments and factorial analyses were designed to isolate and quantify the individual effects of temperature, precipitation, and CO 2 changes, and their interactive effects on NPP.The research objectives are: (1) projecting the trends of urban ecosystem NPP so as to evaluate the performance (growth, coverage, health, etc.) of the urban vegetation in the Phoenix city in the 21st century under different future RCP scenarios; (2) identify the major environmental controls (climate factors or CO 2 ) over the NPP dynamic in the 21st century; and (3) identify the land-use hotspots and sensitive ecosystems in the city, whose NPP may change dramatically in the 21st century.The major purpose is to subsequently understand NPP changes in response to environmental variability and explore urban ecosystem services and urban ecosystem sustainability.

Study Area
The study area (~2902 km 2 made up of 11,608 grids points at 500-m resolution) included Phoenix, Arizona and portions of surrounding agricultural land and desert (central latitude/longitude: 33.52 • /−12.08 • ; average elevation: 340 m) (Figure 1).The city has a population of more than 4 million people (US Census Bureau, 2010).The mean annual precipitation is 180-210 mm, about 65% of which occurs in the winter from November to April and the remaining ~35% occurs as monsoonal thunderstorms from June to August.Mean annual temperature is 22 • C, with hot summers and mild winters.Native vegetation is desert shrub/scrub communities.Urbanized area is occupied by large areas of either cultivated grass and broadleaf trees or desert-like landscaping with drought-tolerant shrub species and gravel ground cover [21] (Table 1).vegetation is desert shrub/scrub communities.Urbanized area is occupied by large areas of either cultivated grass and broadleaf trees or desert-like landscaping with drought-tolerant shrub species and gravel ground cover [21] (Table 1).

Model Description
The HPM-UEM (Hierarchical Patch Mosaic-Urban Ecosystem Model) is a multi-scaled and process-based terrestrial ecosystem model, which explicitly treats spatial pattern and hierarchical structure of urban landscape by incorporating both top-down controls and bottom-up mechanisms in urban environment [21].Based on the conceptual model of Hierarchical Patch Dynamics [22][23][24], the HPM-UEM includes six nested hierarchical levels: region → landscape → land-use → land-cover

Model Description
The HPM-UEM (Hierarchical Patch Mosaic-Urban Ecosystem Model) is a multi-scaled and process-based terrestrial ecosystem model, which explicitly treats spatial pattern and hierarchical structure of urban landscape by incorporating both top-down controls and bottom-up mechanisms in urban environment [21].Based on the conceptual model of Hierarchical Patch Dynamics [22][23][24], the HPM-UEM includes six nested hierarchical levels: region → landscape → land-use → land-cover (or local ecosystem) → plant population → individual plant.The biophysical and ecophysiological processes at and below local ecosystem scale are simulated at daily time-step based on knowledge from natural ecosystems, with carbon and water cycles coupled [25,26].At higher levels, factors associated with landscape patterns, land management, and multiple environmental changes (e.g., climatic conditions, CO 2 ) were explicitly incorporated into the spatially nested land hierarchy.To address the spatial heterogeneity in the landscape and environments, urban lands are treated as spatially nested patch mosaics in the HPM-UEM.By addressing the six hierarchical levels from individual plant to urbanized region, the model provides a "hierarchical ladder" to scale up local ecosystem functions across the nested urban land hierarchies and facilitate linking ecosystem processes and socioeconomic drivers.The model has been extensively evaluated against observed NPP and carbon pools in our study area in a previous study [21].This research was based on the model parameters developed in the previous study [21], with updated land-use map and land-use/land-cover parameters for the Phoenix city, and extended climate and CO 2 data from 2006 to 2099.
In HPM-UEM, NPP (g C/(plant•year)) of individual plant is calculated as where GPP is the gross primary productivity (g C/(plant•year)), and R a is the autotrophic respiration rate (g C/(plant•year)).
The model estimates the photosynthesis rate of a plant (A) following a biochemical model of leaf photosynthesis originally developed by Farquhar et al. (1980) [27] and subsequently expanded by Collatz et al. (1992) [28] and other researchers (Sellers et al., 1996;Bonan, 1996) [29,30].Unlike empirical models, HPM-UEM does not directly simulate the correlation between environmental factors (CO 2 , temperature, precipitation, etc.) on photosynthesis.Instead, photosynthesis is explicitly connected to stomatal conductance (i.e., gs) and controlled by multiple environmental factors.The photosynthesis rate is initially calculated on the leaf level.A = min(w c , w j, w e ) PFT)   where c i and o i are the partial pressure of internal leaf CO 2 and O 2 , respectively (Pa); Pressure is the atmospheric pressure (Pa); Γ * is the CO 2 compensation point (Pa); K c and K o are the Michaelis-Menten constants for CO 2 and O 2 ; absorbed PAR (Rad abs,PAR ; W/m 2 ) is converted to photosynthetic photon flux by assuming 4.6 µmol photons per Joule; α is the quantum efficiency; and V max is the maximum rate of carboxylation varied with temperature (T, • C), and the water potential of the crown (ψ crn , kPa): where V max25 is the maximum carboxylation rate at 25 • C; α vmax is a temperature sensitivity parameter, indicating the magnitude of change in V max with temperature altering every 10 • C away from the reference temperature (25 • C); f(T) is an empirical function that delineates the response of leaf carboxylation to temperature; and f(ψ crn ) is an empirical function that that delineates the response of leaf carboxylation to the leaf water potential.The assimilation rates (GPP; g C/(s m 2 LA)) on a unit-projected leaf area basis for both C 3 and C 4 plants are estimated independently for the sunlit and shaded canopy fractions.The GPP of a plant is then calculated by multiplying A with the leaf area of the plant: where dayl (s) is the length of day and f GPP,SW is the effect of soil water deficiency on productivity.
where tran (mm/s) and Ptran (mm/s) are actual and potential transpiration rate of the plant.The details of the water-uptake and transpiration module has been described by Zhang et al. [21].
where R m (g C/(plant day)) is the daily maintenance respiration calculated as a function of the ambient temperature [21].After maintenance respiration is subtracted from the GPP, 25% of the remainder is taken as growth respiration, which is the cost of producing new tissues [32].

Input Data Description
The model inputs include 500-m resolution spatial datasets of climate, atmospheric CO 2 , land-use and land-cover, and base maps for the study area.The base map datasets contain: (a) the topographic maps (elevation, slope, and aspects) derived from the 7.  [33][34][35].A detailed description of sources and the development of the base map datasets are found in Zhang et al. [21].
Future climate trends were based on the fifth phase of the Coupled Model Intercomparison Project (CMIP5) projections [36] (https://esgf-data.dkrz.de/search/esgf-dkrz/).The CMIP5 strategy includes several types of long-term climate change modeling experiments.The core simulations within the suite of CMIP5 long-term experiments include future projection simulation forced with specified concentrations [referred to as "representative concentration pathways" (RCPs)].The labels for the RCPs provide a rough estimate of the radiative forcing in the year 2100 (relative to preindustrial conditions).For example, the radiative forcing in RCP8.5 increases throughout the twenty-first century before reaching a level of about 8.5 W m −2 at the end of the century.In addition to this "high" scenario, RCP4.5 is intermediate, and there is a low, so-called peak-and-decay, scenario, RCP2.6, in which radiative forcing reaches a maximum near the middle of the twenty-first century before decreasing to an eventual nominal level of 2.6 W m −2 [36,37].The future climate projections from different general circulation model (GCM) models could be very different.Because of the large number of simulation unit (over 11 thousand grid pixels) and long simulation time-period (100 years with a daily time-step) in this study, it was practically impossible to conduct multiple simulations for each RCP scenario (currently, each round of simulation takes more than one month, not including the additional time for data processing).Therefore, different GCM model products were evaluated using historical (2000)(2001)(2002)(2003)(2004)(2005) climate datasets developed based on observational climate records from local meteorological stations [21] (Table 1).Then, based on the evaluation result, the best GCM product was selected to develop the high-resolution future climate dataset to drive the model simulation.In total, there are only eight GCM model products meet the data requirement of this study (Table 1).Our evaluation indicated the MPI-ESM-MR model products have the best performance in the study area, with highest correlation (R 2 = 0.85) and least bias of −81.Therefore, the future climate and CO 2 datasets for the Phoenix (2006-2099) were developed based on the MPI-ESM-MR model products and the RCP protocol, which were statistically downscaled to 500-m resolution based on a high-resolution historical (2000-2005) climate and CO 2 dataset of the Phoenix city [21] using the Localized Constructed Analogs (LOCA) method [38,39].The LOCA spatial downscaling procedure includes four steps: (a) develop analog pool points and spatial masks; (b) select analog days at the regional scale; (c) find the one best matching analog day at the local scale; and (d) construct final downscaled field using scale factors.Another detailed description of the LOCA methodology can be found in Pierce et al. [40].The land-use map used in this study was based on the 2012 Maricopa land-use map developed by the Maricopa Association of Governments (MAG) (http://geo.azmag.gov/maps/landuse).In total, 11 land-use types were identified in the Phoenix city including commercial area, industrial area, mesic residential area, xeric residential area, institutional area, transportation, waters, recreation park, construction transition, agricultural land, and vacant nature, respectively (Figure 1a).To define the land-cover composition of the land-use types, 50 sampling plots (each with 100 m × 100 m size) for each of the 11 land-use types were selected randomly.Then, by visual classification based on a 0.25-m remote sensing images from Google-Earth (Spot5, Phoenix, AZ, USA), we estimated the land fraction of seven land-cover types (woodland, lawn, shrub/bush (shrubland), cropland, waters, bare soil, and impervious surface area) in each land-use sampling plot, and calculated the mean land-cover composition in each of the 11 land-use type in the Phoenix city (Table 2; Figure 1b).

Experimental Design and Factorial Analysis
The implementation of HPM-UEM simulation includes three steps: (a) equilibrium run; (b) spin-up; and (c) transient run.The main purpose of equilibrium run was to establish a baseline for the biomass, soil carbon and water pools.The model was driven by the initial climate datasets and CO 2 concentration in the year of 2000 at the equilibrium step.Then, a spin-up run of 1000 (200 spins × 5 years/spin), driven by random data sequence of de-trend climate data of 2000-2005 was conducted to reduce biases in the transient simulations.Finally, the time series data for climate and CO 2 was applied to simulate the NPP dynamics in the 21st century (2006-2099) under different RCP scenarios.
To identify the magnitude and relative contribution of various environmental factors (temperature, precipitation, and CO 2 concentration), four numeric experiments were designed.The OVERALL experiment (multifactor) was set up to simulate the combined effect of climate and CO 2 change on NPP during the entire study period.The other experiments-CO 2 , TMP, and PPT-were set up to project the relative devotion of CO 2 change (or CO 2 fertilization effect), temperature change, and precipitation change on NPP during the entire period, respectively.In the multifactor experiment (i.e., OVERALL), all environmental factors were allowed to change from day by day.In the three other factorial experiments, one factor was changed while the others were kept unchanged (to the level of equilibrium year); taking the CO 2 experiment for example, CO where E f f ect f actor is the magnitude of NPP change during a certain time period as defined by the subscriptions of "start" to "end" in response to the changes in climate/CO 2 factors (including CO 2 , temperature, precipitation).Following previous studies [41,42], the NPP dynamics during the first-half, second-half, and the whole 21st century were estimated by comparing the five-year averaged NPP of the early-21st century (2006-2010), mid-21st century (2048-2052), and late-21st century (2059-2099).For example, the "start" was set to the 2006-2010 period, and "end" was set to "2048-2052", when estimating the NPP change during the first-half of the 21st century.
In order to compare our simulated CO 2 fertilization effect with other studies, we calculated its growth factor (or beta factor) based on Keeling's formula [43]: where β is the normalized CO 2 effect; and E and B denote experiments under elevated CO 2 and baseline CO 2 , respectively.

Temporal Change Pattern of NPP
During the 21st century (comparing the early-21st century (2006-2010) with the late-21st century (2095-2099)), the regional averaged air temperature in Phoenix city would increase by 0.80 • C, 0.69 • C, and 4.45 • C under RCP2.6,RCP4.5, and RCP8.5 scenarios, respectively; the precipitation would increase 114.3% under RCP8.5 scenario but decrease by 1.9% and 0.2% under RCP2.6 and RCP4.5 scenarios, respectively.The atmospheric CO 2 concentration increased by 40.5 ppm, 156.99 ppm, and 545.85 ppm under RCP2.6,RCP4.5, and RCP8.5 scenarios, respectively, as shown in Figure 2. The simulation indicated the climate changes and elevated CO 2 concentration would result in a net increase of urban ecosystem NPP in the Phoenix city (from 2006-2099), as shown in Table 3 and Figure 3a.The increase rates of NPP would differ strongly under the different future climate scenarios, ranging from 15% under the RCP2.6 scenario to 100% under the RCP8.5 scenario.During the 21st century, the NPP would increase steadily under both the RCP4.5 and RCP8.5 scenarios, although the NPP increase rate of the later would be twice higher than the former.Under the RCP8.5 scenario, the NPP would increase 41% in both the 1st half and the 2nd half of the 21st century.Under the RCP4.5 scenario, the increase rate of NPP would slow down slightly from 25% during the 1st half of the century to 21% during the second half of the century, possibly related to the declined precipitation.In contrast, the NPP would first increase 16%, reaching peak by the mid of the 2040s, after which the growing trend would disappear under RCP2.6,possible related to the declined precipitation and the atmospheric CO 2 concentration.The NPP changes in the 1st half of the 21st century was calculated as the difference between the mean NPP in the mid-21st century and that in the early-21st century.The NPP changes in the 2nd half of the 21st century was calculated as the difference between the mean NPP in the late-21st century and that in the mid-21st century.The NPP changes in the 21st century was calculated as the difference between the mean NPP in the late-21st century and that in the early-21st century.Unit: K ton C (1 K = 1000).

Relative Contributions of Environmental Factors to NPP Changes
The results showed that relative contribution of environmental factors on urban ecosystem NPP would vary significantly in Phoenix (Table 4).The overall trends of the urban NPP would apparently be dominated by the CO2 fertilization effects, which would account for 56-61% of the NPP increase during the 21st century (Table 4).The dominant role of the CO2 effect on future NPP would be

Relative Contributions of Environmental Factors to NPP Changes
The results showed that relative contribution of environmental factors on urban ecosystem NPP would vary significantly in Phoenix (Table 4).The overall trends of the urban NPP would apparently be dominated by the CO 2 fertilization effects, which would account for 56-61% of the NPP increase during the 21st century (Table 4).The dominant role of the CO 2 effect on future NPP would be particularly obvious under the RCP4.5 and RCP8.5 scenarios (Figure 3c,d).The air temperature and precipitation would influence the inter-annual fluctuation of the NPP (Figure 3b-d).The temperature change would contribute to 49% of the NPP increase under the RCP2.6 scenario (Table 4).The precipitation change would contribute to 20% of the NPP increase under the RCP4.5 scenario (Table 4).It is noteworthy that the reduced precipitation (19%) would result in a 4% loss in NPP during the first half of the 21st century, and the declined CO 2 concentration would lead to a 4% loss in NPP during the 2nd half of the 21st century, under the RCP2.6 scenario (Table 4; Figure 3e).

Spatial Variation in NPP and the Responses from Different Ecosystem Types
Under the RCP2.6 scenario, the NPP of about 39.38% of the urban ecosystem, mostly located in the western and southern Phoenix, would not change much (the grey areas in Figure 4a) during the 21st century, while the strongest NPP increase (generally >60 g C/m 2 ) would be found in the agricultural lands located in western and southeastern Phoenix.Under the RCP4.5 and RCP8.5 scenarios, most of the urban ecosystems, except for the industrial and transportation lands that mainly locate to the north of Salado River, would see a significant increase of NPP (>40 g C/m 2 under the RCP4.5 scenario, and >75 g C/m 2 under the RCP8.5 scenario) during the 21st century, with the strongest NPP increase (>200 g C/m 2 under the RCP8.5 scenario) found in the mesic residential areas that mainly locate in the east, south, and southwest of the Phoenix Mountains Preserve (including Moon Hill, Shaw Butte, North Mountain, Stoney Mountain, Dreamy Draw, Squaw Peak, and the informally named Quartzite Ridge [44]) (Figure 4b,c).
Among the different ecosystem types (or land-cover type), the crop ecosystems would have the strongest NPP increase (at least 43% stronger than any other ecosystem type), followed by the urban lawn and forests, while the remnant dryland including shrubland and grassland would have the weakest NPP increase, under the RCP4.5 and RCP8.5 scenarios (Figure 5).Under the RCP2.6 scenario, the situation will be more complicated.The NPPs of all ecosystem types were projected to increase during the 1st half of the 21st century.Then, as the CO 2 concentration declining in the 2nd half of the century under the RCP2.6 scenario, the increase trend of NPP would disappear (in the natural dryland ecosystems) or even reverse (in the urban forest and lawn that would have accumulated large biomass during the first half of the century).

Discussion
During the 21st century, the urban ecosystem NPP would universally increase in Phoenix city, whereas the magnitude of NPP trend was different under RCP2.6,RCP4.5, and RCP8.5 scenarios.It is noteworthy that the NPP growth rate exhibited significant variability (i.e., NPPRCP8.5 > NPPRCP4.5 > NPPRCP2.6) with an increase in radiative forcing, as shown in Figure 4. Recent research shows that different radiative forcings have markedly contrasting physiological effects on terrestrial ecosystems [45,46].Another land surface model, JULE (The Joint UK Land Environment Simulator), has shown that a rise in atmospheric CO 2 concentration ([CO 2 ]) from 354 to 426 ppm corresponds to an increase in radiative forcing 2 , which generated a higher NPP trajectory [47].Our simulation also indicates that the effect of CO 2 fertilization on NPP outcompetes that of other environmental factors (such as temperature and precipitation) in the 21st century.Zhao et al. found that higher CO 2 concentration would strongly enhance terrestrial NPP during the growing season, based on 10 Earth system models participating in the CMIP5 Project (i.e., fifth phase of the Coupled Model Intercomparison) [48].This is broadly consistent with experimental evidence of CO 2 fertilization of photosynthesis [49] and an increase in water use efficiency under elevated CO 2 concentration [50].While the HPM-UEM-simulated β-factor (i.e., the normalized CO 2 fertilization effects) fell within the ranges of field observations for grassland and forest, the model underestimated the CO 2 fertilization effect on shrubland/desert by 40% to 53%, similar to another model study [51] (Table 5).In addition, our study predicted a 2-4% increase in crop NPP in response to each 100 ppmv [CO 2 ] increase, matching well with the estimates (2.4-4.2%) of Fischer (2009) [52] but lower than the estimates (5-17%) of McGrath and Lobell (2013) [53].Therefore, the model simulations might have underestimated the impacts of CO 2 elevation on the NPP of the remnant shrubland/desert in the Phoenix city.As a key indicator of ecosystem function, NPP has been shown to correlate with the overall value of ecosystem services [67].According to Costanza et al. (2007), a 1% increase in NPP corresponding to 2.6% increase in ecosystem value [68].The significant increases of the NPP under the RCP4.5 and RCP8.5 scenarios imply the climate and CO 2 changes may enhance the urban ecosystem services in the 21st century.In particular, the simulated results indicate the cropland yield may increase rapidly in the 21st century (Figure 4).Other studies also showed future climate change may benefit food provision [69,70].Improved plant growth and leaf area will also enhance other urban ecosystem services such as heat regulation [71], buffering capacity or the absorption capacity for wastes and emissions [72,73], and carbon sequestration [74].These ecosystem services would partially moderate the worsen heat stress and air pollution in urban areas under the future scenarios [75][76][77][78][79].It is noteworthy that the wealthy population living in the mesic residential areas may enjoy more benefit from the enhanced ecosystem services than the poorer population living in the xeric residential areas, because the NPP of the mesic ecosystems (lawns and urban forests) would increase much more quickly than the xeric ecosystems (dry shrubs and grassland) (Figures 4 and 5).
Although our model projections indicate that higher atmospheric CO 2 may stimulate ecosystem NPP in the urban area, it does not mean that a high CO 2 emission scenario is preferable, nor does it imply that the enhanced ecosystem productivity could offset the elevated urban CO 2 emission.Even if not taking into account the higher soil CO 2 emission due to enhanced heterotrophic respiration under a warmer climate [80], the magnitude of ecosystem NPP in the Phoenix city (629 k ton C/year by the late-21st century, under the RCP8.5 scenario) is too small compared with the anthropogenic CO 2 emission (71,356 K ton C/year according to Koerner and Klopatek, 2002) [81].In fact, the 20 largest U.S. cities each year contribute more CO 2 to the global atmosphere than the total land area of the continental United States can absorb [82].
It also should be aware that NPP is only one of the many indicators of ecosystem health.Increased NPP would not necessarily means improved ecosystem sustainability.Rapid climate changes may benefit some plant species more than others (Figure 4) and trigger ecosystem succession [83][84][85].Since the desert cities like Phoenix are hotspot of climate changes, it is possible that some native species will be at a disadvantage in face of the competition from alien species which can take full advantage of the elevated CO 2 in the future [86][87][88][89][90].The impacts of future climate change on the ecosystem diversity and native ecosystem conservation in the Phoenix metropolitan area are still unclear.
The enhanced vegetation growth in Phoenix may partially alleviate the increased heat stress and pollution under future climate change scenarios.Because of intensive human managements, urban ecosystems may perform better than the natural ecosystems in response to climate stresses [91].The enhanced NPP of the desert urban ecosystem will cost large amount of irrigation water.As more water resources will be used in urban irrigation, the natural desert ecosystems around the city might suffer water shortage in the future [92][93][94].In other words, we may expect an improved ecosystem NPP in Phoenix at the cost of serious natural ecosystem degradation in this arid region in the future.
Our study may also underestimate the negative effects of future climate change on ecosystem NPP.A recent study found that climate extremes may reduce the CO 2 fertilization effect [95].Because the HPM-UEM model, like most general ecosystem models, was designed and parameterized to simulate the ecological processes under ordinary climate condition, its performance under climate extremes are questionable.It is possible that the model may underestimated the negative effects of future climate change on NPP.Moreover, the model does not consider the effects of climate change on tropospheric ozone concentration, diseases and disturbances (e.g., increasing wild-fire events), which may threaten the ecosystem sustainability under the future scenarios [53,[95][96][97][98][99][100][101]].An ecosystem with higher NPP is not necessarily more resistant to the increased disturbances (drought and flood) in the future.Higher vegetation density plus higher temperature and moisture in urban area may stimulate disperse of pests and infectious diseases [102].Our study indicates that the rich (those live in mesic residential area) would enjoy more benefits from the improved urban ecosystem services (e.g., better tree shading and a greener landscape) than the poor (those live in xeric residential area and work in open areas such as transportation, construction areas, etc.).However, the underprivileged population is more vulnerable to pests and infectious diseases.Therefore, it is necessary to allocate more resources to the underprivileged population to help them cope with the heat stress (will rise 4.5 • C by the end of the century) and threats from pests and infectious diseases due to the rapid climate change.
Finally, because of the limitation in time and resources (see Section 2.3), this study only used one model driven by one suite of future climate dataset.Although the climate dataset matches well with the historical climate records (Table 1), it does not necessarily mean it can predict future climate change accurately.Multiple ecosystem models with multiple climate datasets approach could reduce the uncertainties in future projections (Wallach et al., 2016) [103].Meanwhile, uncertainties may exist in our estimated NPP dynamics, which were calculated by comparing the five-year mean NPPs of the beginning and ending years of the study periods.Although this is a common practice in previous studies [41,42], climate anomalies in these years may lead to miscalculations of NPP change trends.

Conclusions
Urban vegetation provides ecological services that promote both the ecosystem integrity and human well-being of urban area, and thus is critical to urban sustainability.Using NPP as an indicator, this study investigated the impacts of future climate change on the health of urban vegetation and urban sustainability.A process-based urban ecosystem model, HPM-UEM, was introduced to explore spatiotemporal dynamic of urban ecosystem NPP in the Phoenix city (AZ, USA) under RCP2.6,RCP4.5, and RCP8.5 scenarios during the 21st century.The simulated results found that urban ecosystem's NPP would increase from 2006 to 2099, whereas NPP variations exhibited obvious differences in the amounts and spatial distributions under different RCP scenarios.The magnitude of NPP increase by 14% (RCP2.6),51% (RCP4.5),and 99% (RCP8.5)compared to that in the late 2000s.According to factorial analysis, CO 2 fertilization effect on NPP variation (accounting for 56-61%) outcompetes that of climate change in the 21st century.NPP trend in Agricultural land and mesic residential area would increase considerably under three scenarios, primarily owing to elevated CO 2 concentration and increase in water use efficiency.Although higher ecosystem NPP in the future implies improved ecosystem services that may help to alleviate the heat stress (by providing more shading) and air pollution in the city, this will be at the cost of higher irrigation water usage, probably leading to water shortage in the natural ecosystems in this arid region.Our study indicated the rich (such as in mesic residential area) would enjoy more benefit from the improved urban ecosystem services than the poor (such as in xeric residential area).

Figure 1 .
Figure 1.Hierarchical land structure of Phoenix, AZ USA: (a) the study area of Phoenix and land use composition of the urban landscape; and (b) land cover composition of land use functional types.

Figure 1 .
Figure 1.Hierarchical land structure of Phoenix, AZ USA: (a) the study area of Phoenix and land use composition of the urban landscape; and (b) land cover composition of land use functional types.
5-min US Geological Survey (USGS) National Elevation Dataset (NED); (b) the soil maps contain bulk density (g/cm 3 ), volumetric content of sand and clay (%), depth (m), and acidity (pH), which are generated from 1-km resolution digital general soil association map; and (c) spatial map of nitrogen deposition based on the study of Grossman-Clarke et al. (2005), Baker et al. (2001), and Lohse et al. (2008) 2 concentration was transient during 2000-2099 while climate data (temperature and precipitation) were kept at the equilibrium level.Based on the simulation results, a serial of factorial experiments were conducted to quantify the individual effects of CO 2 /climate factors: E f f ect f actor = NPP f actor, end − NPP f actor, start (where factor [CO2, TMP, PPT])

Figure 3 .
Figure 3. Influence of environmental factors on NPP change in Phoenix since early 21st century as simulated by the HPM-UEM under RCP2.6,RCP4.5, and RCP8.5 scenario, respectively: (a) A fiveyear moving average NPP of multiple factors (temperature, precipitation and CO2) under RCP2.6,RCP4.5, and RCP8.5 scenario, respectively.(b-d) Changes in NPP since the early 21st century (i.e., in comparison to the mean annual NPP of 2006-2010) in response to different factors under RCP2.6,RCP4.5 and RCP8.5 scenario, as predicted by the numeric experiments: CO2 is the carbon dioxide change only experiment; TMP is the temperature change only experiment; PPT is the precipitation change experiment; and the OVERALL experiment simulate the synergistic effects of CO2, temperature and precipitation.

Figure 3 .
Figure 3. Influence of environmental factors on NPP change in Phoenix since early 21st century as simulated by the HPM-UEM under RCP2.6,RCP4.5, and RCP8.5 scenario, respectively: (a) A five-year moving average NPP of multiple factors (temperature, precipitation and CO 2 ) under RCP2.6,RCP4.5, and RCP8.5 scenario, respectively.(b-d) Changes in NPP since the early 21st century (i.e., in comparison to the mean annual NPP of 2006-2010) in response to different factors under RCP2.6,RCP4.5 and RCP8.5 scenario, as predicted by the numeric experiments: CO 2 is the carbon dioxide change only experiment; TMP is the temperature change only experiment; PPT is the precipitation change experiment; and the OVERALL experiment simulate the synergistic effects of CO 2 , temperature and precipitation.

Figure 4 .Figure 4 . 20 Figure 4 .Figure 5 .
Figure 4. Spatiotemporal pattern of NPP (g C/m 2 ) in the Phoenix city during 2006-2099 as simulated by HPM-UEM: (a) RCP2.6 scenario; (b) RCP4.5 scenario; and (c) RCP8.5 scenario.NPP change of 2006-2099 is the value of late 21st century subtract that of early 21st century; NPP of early 21st century is the mean value during 2006-2010; NPP of late 21st century is the mean value during 2095-2099; NPP is net primary productivity; RCP is representative concentration pathway; 2.6, 4.5 and 8.5 are the radiative forcings of year 2100 relative to that of year 1750.
max 0, (ψ crn − ψ close ) ψ open − ψ close } where T min , T max and T opt are minimum, maximum, and optimum temperature ( • C) for photosynthesis, respectively; ψ close and ψ open are thresholds of crown water potential (kPa) at which stomata of leaf begin to fully close and fully open, respectively; and ψ crn is calculated as

Table 1 .
Evaluating the CMIP5 GCM products using historical (2000-2005) monthly climate dataset developed based on field observations by Zhang et al. [21].

Table 2 .
Description of the 11 land-use types in the Phoenix city.

Table 3 .
Model projected changes in the total ecosystem NPP of the Phoenix city, AZ, during the 21st century.
Note: RCP is Representative Concentration Pathway, and 2.6, 4.5 and 8.5 are radiative forcing of year 2100 in contrast to year 1750.Mean annual NPP of the early-, mid-and late-21st century are the mean NPPs of 2006-2010, 2048-2052 and 2095-2099, respectively.

Table 4 .
Results of factorial analysis in Phoenix.: the multiple factors experiment that simulated combined effect of climate and CO 2 change on NPP; CO 2 : the CO 2 experiment that CO 2 was transient during 2000-2099 while climate data were equilibrium level; TMP: the temperature experiment that temperature was transient while CO 2 and precipitation were equilibrium level; PPT: the precipitation experiment that precipitation was transient while CO 2 and temperature were equilibrium level; and INTACTIVE: the interactive effect between climate change and CO 2 change.RCP2.6:RCP is Representative Concentration Pathway, while 2.6, 4.5 and 8.5 are the radiative forcings of year 2100 in contrast to year 1750.NPP is net primary productivity.Unit: K ton C. OVERALL

Table 5 .
Comparison of NPP trend derived from CO 2 change between HPM-UEM model and FACE/OTC experiment (FACE: Free-air CO 2 enrichment).