Precipitation response to climate change and urban development over the continental United States

Appropriately characterizing future changes in regional-scale precipitation requires assessment of the interactive effect owing to greenhouse gas-induced climate change and the physical growth of the built environment. Here we use a suite of medium resolution (20 km grid spacing) decadal scale simulations conducted with the Weather Research and Forecasting model coupled to an urban canopy parameterization to examine the interplay between end-of-century long-lived greenhouse gas (LLGHG) forcing and urban expansion on continental US (CONUS) precipitation. Our results show that projected changes in extreme precipitation are at least one order of magnitude greater than projected changes in mean precipitation; this finding is geographically consistent over the seven CONUS National Climate Assessment (NCA) regions and between the pair of dynamically downscaled global climate model (GCM) forcings. We show that dynamical downscaling of the Geophysical Fluid Dynamics Laboratory GCM leads to projected end-of-century changes in extreme precipitation that are consistently greater compared to dynamical downscaling of the Community Earth System Model GCM for all regions except the Southeast NCA region. Our results demonstrate that the physical growth of the built environment can either enhance or suppress extreme precipitation across CONUS metropolitan regions. Incorporation of LLGHGs indicates compensating effects between urban environments and greenhouse gases, shifting the probability spectrum toward broad enhancement of extreme precipitation across future CONUS metropolitan areas. Our results emphasize the need for development of management policies that address flooding challenges exacerbated by the twin forcing agents of urban- and greenhouse gas-induced climate change.


Introduction
Appropriately characterizing trends in future precipitation has implications for ecosystems and societies across the world. For example, plants and animals accustomed to an envelope of climatic variability must maintain poleward movement along the planet's surface to ensure pace with historically experienced environmental conditions (Loarie et al 2009). Likewise, as populations grow, rapid groundwater extraction (Scanlon et al 2012, Dalin et al 2017 and freshwater scarcity (Mekonnen and Hoekstra 2016) place undue pressure on billions of people across the world. While management of water resources depends on efficient use, infrastructure investment and sound policies, characterizing changes in future mean and extreme precipitation is similarly important.
Global climate model (GCM) simulations indicate that the poleward expansion of the subtropical dry zones resulting from increased emissions of long-lived greenhouse gases (LLGHGs) is a key driver of future precipitation change (Seager et al 2007), though uncertainties about the rate of expansion remain (Quan 2014). Increased atmospheric concentration of LLGHGs may also drive future changes in precipitation through multiple land-ocean-atmosphere processes working in concert. For example, the release of excess latent heat, due to stronger western tropical Pacific storms whose origin is tied to increased sea surface temperatures, has been shown to interact with and amplify midlatitude planetary waves. Coupled with a weaker jet stream, this may lead to increased persistence of storm systems and extreme precipitation (Francis and Vavrus 2012, Muller 2013, Palmer 2014, Guilbert et al 2015. Appropriately projecting effects of extreme precipitation is particularly important because of widespread societal implications including flash flooding, landslides in complex terrain, degradation of aquatic ecosystems resulting from enhanced stormwater runoff, and extensive damage to infrastructure and loss of life. For example, some locations along Colorado's Front Range experienced as much rainfall during the 2013 September floods as normally received during an entire year (Gochis et al 2015). The stalled cold front and ensuing extreme rainfall became the state's costliest natural disaster, resulting in upwards of $1 billion in property damage (Aguilar 2017).
There is mounting evidence that urban areas also modify precipitation (Changnon et al 1971, Changnon 1979, Lowry 1998, Bornstein and Lin 2000, Rosenfeld 2000, Shepherd 2005, Liu and Niyogi 2019. Destabilization of the urban boundary layer resulting from development of the urban heat island (UHI) circulation, modification of convective structures due to the rough urban form, and the presence of cloud condensation nuclei in more polluted urban environments are broadly considered responsible for urban-induced precipitation modification. Satellite and radar data have provided additional means of identifying urban-induced precipitation anomalies (e.g. Rosenfeld 2000, Shepherd 2006, Mote et al 2007, Kishtawal et al 2010, Ashley et al 2012, Haberlie et al 2015Wu et al 2019. Ensuring that identification of urban-induced precipitation anomalies are physically based requires numerical approaches (i.e. processbased models) to confirm the role of the urban environment in altering precipitation amount, pattern and timing (Pielke et al 2007, Georgescu et al 2009, Li et al 2013, Benson-Lira et al 2016, Paul et al 2018, Zhang et al 2018, Kusaka et al 2019, Yang et al 2019, Zhang et al 2019. Use of regional climate models (RCMs), due to their higher spatial and temporal resolution, enables more detailed process-based understanding and simulation of climate at regional scales (Rummukainen 2010).
The processes driving precipitation are complex and occur at multiple spatial and temporal scales. Physical mechanisms associated with increased concentration of LLGHGs are fundamentally different than those associated with the built environment. For example, Clausius-Clapeyron dependent scaling relationships that relate temperature to extreme precipitation are often used to estimate the influence of a warming world on precipitation extremes (Trenberth et al 2003). However, such scaling relationships assume a constant relative humidity, no changes in the general circulation and indicate super (i.e. greater than) or sub-scaling depending on distinct thermal regimes (Panthou et al 2014, Westra et al 2014, Ivancic and Shaw 2016, Giorgi et al 2019. Development of average and extreme future precipitation projections at regional scales therefore requires explicit accounting of large-scale atmospheric dynamics and thermodynamics associated with increasing greenhouse gas concentrations. In addition, precipitation projections must incorporate effects arising from the rapidly expanding built environment. In particular, accounting for processes associated with local-scale thermodynamic/dynamic impacts caused by the built environment and interaction with LLGHG-driven thermodynamic changes and related modification of large-scale circulation is essential. While the urbaninduced thermal effect has been shown to be of similar order of magnitude as that owing to twenty-first century greenhouse gas-induced climate change on regional urban scales (Georgescu et al 2013, 2014, Krayenhoff et al 2018, the interplay between these forcing agents remains unknown with respect to precipitation. The interactive effect between greenhouse gasinduced climate change and the growth of the built environment on regional-scale precipitation has received scant attention to date. To our knowledge there is one exception, although it assumes changes in large-scale circulation patterns due to global warming are independent from changes arising from urbanization (Shastri et al 2019). Here we use a suite of decadal scale RCM simulations that include a contemporary and a range of potential future continental US (CONUS) climate scenarios with varied LLGHG forcing and urban expansion scenarios. We dynamically downscale two GCMs and their respective Representative Concentration Pathways (RCPs) with the Weather Research and Forecasting (WRF) model to build on previous research in two key ways: (a) demonstrate WRF skill in simulating the spatial pattern and magnitude of observed mean and extreme seasonal evolution of CONUS precipitation; (b) examine projected seasonal extreme precipitation change across CONUS and the corresponding contribution from the twin forcing agents of urbanization and LLGHGs.
This work is, to our knowledge, the first attempt to generate regional-scale projections of precipitation change through the explicit assumption that the twin forcing agents of urbanization and LLGHGs interact with one another. Although our focus is over CONUS, this work should motivate similar analyses across regions of the world where these agents are expected to play an increasingly important role in modifying regional climate, including, for example, the rapidly growing urban regions of Asia and Africa (Jones and O'Neill 2016).

Experimental design
We utilize decadal scale WRF RCM (version 3.6) simulations focused over CONUS to examine the response of mean and extreme precipitation to endof-century climate change and urban expansion using a suite of time-slice mode numerical experiments (see Krayenhoff et al 2018). A continuous WRF simulation representing contemporary climate-hereafter Control-is conducted for 2000-2009. The Control experiment is initialized and forced at the lateral boundaries using data obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-interim reanalysis (ECMWF 2009). WRF is coupled to a single-layer urban canopy model (SLUCM) which represents fundamental physical processes associated with the built environment and its interaction with the overlying atmosphere (Kusaka et al 2001, Chen et al 2011. To account for increased concentration of LLGHGs, WRF simulations dynamically downscale RCP 8.5 from two GCMs for the end of the current century: 2090-2099. Specifically, we use (a) the biascorrected Community Earth System Model (CESM) (Monaghan et al 2014) and (b) the Geophysical Fluid Dynamics Laboratory (GFDL) (Dunne et al 2012) model. We select these models, among the suite of available CMIP5 models, because of their varying sensitivities to enhanced greenhouse gas forcing. During each 10 year simulation period, urban extent and density are included based on projections from the Integrated Climate and Land Use Scenarios dataset (Bierwagen et al 2010; see section 2.2 below). Our experimental design permits analysis of precipitation response to increased emissions of LLGHGs resultant from realistic physical expansion of the built environment, both separately and in tandem (table 1). Further details on simulation setup and domain configuration are available in Krayenhoff et al (2018).
The physics options utilized in the WRF simulations described above account for microphysics, long-and shortwave radiation, convection, as well as boundary and surface layer processes. We utilized the unified Noah land surface model (LSM) to represent water and energy exchange between the soil column and overlying atmosphere (Tewari et al 2004). The Noah LSM was utilized over vegetated nonurban grid cells as well as for the pervious portion of urban grid cells. All physics options are detailed in Krayenhoff et al (2018;supp. table 1). In addition, the SLUCM physical parameters (e.g. urban land use class dependent biophysical parameters) are incorporated through WRF input parameter tables. These are also detailed in Krayenhoff et al (2018;supp. table 2). All WRF RCM CONUS simulations are conducted at 20 km grid spacing, with output variables saved at a 3 h frequency. In all, here we examine a suite of six WRF RCM decadal experiments, including the Control simulation and alternative futures accounting for climate change and urban expansion. Collectively, these numerical experiments span more than half a century of simulation time (see Krayenhoff et al 2018;supp. table 3). Given the breadth of decadal scale modeling scenarios conducted, direct dynamical downscaling of more than two GCMs exceeded our computational and storage capability. Nonetheless, in addition to addressing the core scientific questions posed, our results are also expected to highlight simulated uncertainty owing to the use of two independent GCMs.

Datasets used and model evaluation
We incorporate the physical growth of the urban environment using the Integrated Climate and Land Use Scenarios (version 1.3.2; hereafter ICLUS) A2 population projections (Bierwagen et al 2010). We utilize ICLUS to derive spatially explicit impervious fraction and urban density, as an input to WRF and the SLUCM, across CONUS for Control and future dynamically downscaled experiments. Use of ICLUS data permits characterization of end-ofcentury urban development. We incorporate regionally suitable background vegetation for each urban grid cell using Moderate Resolution Imaging Spectroradiometer data; because we do not have appropriate characterization of background vegetation for the end of century, we assume it remains constant in time.
Finally, we utilize the Climate Prediction Center's (CPC's) gridded Unified Gauge-Based Analysis of daily precipitation dataset to examine precipitation skill of the Control simulation. This dataset has a 0.25 • × 0.25 • spatial resolution and was provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA, from their Web site at https:// psl.noaa.gov/data/gridded/data.unified.daily.conus. html. To conduct grid cell by grid cell comparisons, the Control WRF simulation is resampled to the resolution of CPC gridded data using bilinear interpolation. Regions outside CONUS are masked out for purposes of model evaluation and subsequent analyses.

Characterizing the urban-induced signal on extreme precipitation
To improve understanding of urban development and LLGHGs on simulated extreme precipitation, we examine changes in the UHI intensity resulting from expansion of the built environment conditional on changes in daily extreme precipitation. We define extreme precipitation as the 90th, 95th, and 99th percentile of daily mean metropolitan region Control precipitation. The basis for this approach relies on warming-induced destabilization of the urban boundary layer resulting from development of the UHI circulation.
We ascertain the impacts of urban development and LLGHG-induced warming on precipitation at the metropolitan area level by investigating effects across 47 metropolitan regions as defined by Broadbent et al (2020). The metropolitan regions are defined by bounding boxes which enclose end-of-century urban extent of each urban region. The bounding boxes of these regions therefore include a buffer zone surrounding each urban area (see SI appendix, figure S1 in Broadbent et al (2020)) of 20-60 km, which is appropriate for assessing urban-induced precipitation impacts as occurring over, downwind or upwind of the built environment (e.g. Liu and Niyogi 2019).
To determine the change in precipitation, we calculate the 90th, 95th, and 99th percentile of daily mean metropolitan area precipitation (i.e. the daily mean of all pixels in the Control simulation metropolitan area bounding box). On those days corresponding to the 9Xth percentile of daily precipitation, we compute the change in UHI intensity (i.e. ∆UHI) using the following: where T city represents the metropolitan region mean of subgrid urban air temperature and T rural denotes the metropolitan region mean of subgrid rural air temperature for simulation 1 and simulation 2 (i.e. SIM1 and SIM2, respectively). We define subgrid temperature as arising from the separate pervious and impervious portions of the grid square (Krayenhoff et al 2018

Results and discussion
The WRF model has been applied across many regions globally and its precipitation simulation skill has been well documented (e.g. Wang et al 2017, Yang et al 2019. Nonetheless, given the specifics of our domain configuration, a necessary step prior to examining precipitation effects due to increased emissions of LLGHGs and urban expansion is a demonstration of WRF performance for the contemporary (i.e. Control; SIM1 in table 1) climate.

Control precipitation evaluation
Our evaluation of the Control simulation begins with an assessment of averaged and extreme precipitation.
In both instances we characterize WRF performance as annual and seasonal averages. Figure 1 shows the WRF-simulated Control and observed annually averaged mean precipitation and their difference for 2000-2009. WRF captures the principal spatial features of precipitation maxima and minima. The precipitation maxima along the Ohio River Valley, the sharp west-east gradient across the 100th meridian, the reduced precipitation amount across the Intermountain West, and the relatively higher precipitation accumulation across the higher terrain of the Cascades and Sierra Nevada are well represented. Annually averaged differences further highlight the extent of spatial correspondence between the Control and observations, although a broad dry bias of 1-3 mm d −1 is evident along the Gulf Coast states and southern Great Plains states (figure 1(c)). We note that this dry bias is a common feature among RCM simulations (Mearns et   and observed (green bars) distributions. The regional probability distributions highlight distinct precipitation regimes. For example, compare the low peak probability of occurrence and long-tail for the Northwest region to the substantially higher peak probability of occurrence and short tail for the Northeast region. WRF represents these distinct features (e.g. greater frequency of lower precipitation intensity along the Northwest NCA region), demonstrating the model can simulate the precipitation processes and appropriate intensity unique to regional geographies. Notably, WRF demonstrates similar seasonally averaged skill both in terms of spatial representation and shape/magnitude of the probability distribution (supp. figures 1-4 (available online at stacks.iop.org/ERL/16/044001/mmedia)). Figure 2 shows the WRF-simulated Control and observed 99th percentile of daily averaged precipitation and their difference for 2000-2009. The spatial variability of extreme precipitation is captured across the entirety of CONUS. The relative maxima along the eastern seaboard and south-central US is also well represented. Spatial differences between simulated and observed extreme precipitation show but a few regions of bias that are generally small relative to the magnitude of observed and simulated extreme precipitation; for example, a wet bias of 10% or less over the Sierra Nevada and a dry bias of 10% or less over locations along the immediate Gulf Coast areas are apparent (figure 2(c)). Also presented in figure 2 are box and whisker plots illustrating the variability of extreme precipitation (right hand side of spatial panels) for both simulated (pink) and observed (green) distributions for each NCA region. The Control simulation reasonably reproduces the observed variability of extreme precipitation across NCA regions. The large variability in extreme precipitation over the Southwestern NCA region is captured by the Control simulation. Similarly, the observed relatively narrow distribution evident for the Great Plains North and Northeast NCA regions is also captured by the Control simulation. As shown previously for mean precipitation, WRF correctly represents the spatial structure of extreme precipitation across all seasons in terms of both magnitude and geographical extent. Moreover, the distribution of extreme precipitation events is well simulated across NCA regions and seasons (supp. figures 5-8 (available online at stacks.iop.org/ERL/16/044001/mmedia)). Finally, we utilize Taylor diagrams (Taylor 2001) to provide a succinct statistical summary of simulated skill with respect to observed precipitation. Taylor diagrams for annually averaged means and extremes (supp. figure 9 (available online at stacks.iop.org/ERL/16/044001/mmedia)) indicate varying Control simulation skill depending on NCA region. For example, the Southwest NCA region illustrates a correlation coefficient in excess of 0.9, a normalized standard deviation of less than 1.25, and a centered root-mean-squared error of about 0.5 for annually averaged Control precipitation relative to CPC gridded data. As expected, skill is reduced for extreme precipitation. The performance is higher for annual means relative to extremes, and is worst for the convective (i.e. summer) season (supp. figure 10 (available online at stacks.iop.org/ERL/16/044001/mmedia)).
Overall, the Control simulation compares well with observationally based gridded product data and provides confidence in the model's ability to reproduce the seasonally evolving mean and extreme characteristics of precipitation in all CONUS regions. Our results demonstrate mean and extreme precipitation skill that is on par with or better than similar efforts (e.g. Singh et al 2013, Wang and Kotamarthi 2015). The Control simulation provides confidence in WRF's ability to project changes in these metrics for a suite of future scenarios that account for climate change and urban development.

Effects of LLGHG emission and urbanization
We next compare the WRF-simulated change in annual mean and extreme precipitation for each NCA region due to emissions of LLGHGs and urban development. We define the projected change in extreme precipitation as the simulated increase in the 99th percentile of daily averaged Control precipitation. First, we present differences between dynamically downscaled WRF simulations driven by CESM RCP 8.5 and GFDL RCP 8.5 for the end of century and Control (figure 3). Each point in the figure represents a pixel value within the corresponding NCA region. The x-axis denotes the difference of daily averaged precipitation for the 10 year period (i.e. end-of-century decade minus start-ofcentury decade) while the y-axis denotes the difference of 99th percentile of daily precipitation for the 10 year period (i.e. end-of-century decade minus start-of-century decade). Three main conclusions can be drawn: (a) future precipitation changes are projected to be nearly entirely positive, indicating anticipated increases in mean and extreme precipitation across all NCA regions and for both dynamically downscaled GCMs; (b) projected changes in extreme precipitation are greater than projected changes in mean precipitation for all NCA regions (note that red and green lines have slopes that are greater than the slope of the gray line in figure panels, which represents the 1-to-1 correspondence between mean and extreme changes); (c) the ratio of projected change in extreme to mean precipitation is greater for the GFDLrelative to CESM-dynamically downscaled WRF simulation for all regions except the Southeast NCA region, which depicts similar results for both dynamically downscaled GCMs.
Simulated increases in extreme precipitation are at least one order of magnitude greater than projected mean changes, with consistency among NCA regions and the pair of dynamically downscaled GCMs (the lone exception is the CESM RCP 8.5 simulation for the Northeast NCA region; figure 3(d)). These results are broadly consistent across seasons and NCA regions, although exceptions are evident; Table 2. WRF-simulated median change in 99th percentile of precipitation between dynamically downscaled WRF simulations driven by GFDL RCP 8.5 (2090-2099 with urban expansion and Control (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), and CESM RCP 8.5 (2090-2099) with urban expansion and Control (2000-2009) for each NCA region. Units are mm d −1 . Asterisk indicates that the median of simulated change is statistically significantly positive based on one-tailed one-sample Wilcoxon test. Pound denotes the same as the asterisk, but for statistically significantly negative changes. All tests are controlled with familywise type-I error rate less than 0.05 (p-value < 0.0003 for each test under the Bonferroni correction for multiple hypothesis tests). Median changes with opposite signs between CESMRCP85-Urban2100 and GFDLRCP85-Urban2100 are highlighted in bold. However, even in these instances, changes in extreme precipitation are substantially greater than mean changes. Consequently, our focus from this point forward will be on examination of impacts on projected changes in extremes. Seasonally and temporally averaged differences in extreme precipitation for dynamically downscaled simulations driven by end-of-century GFDL and CESM RCP 8.5, both with urban expansion, and Control illustrate broad geographical consistency (figure 4). Both dynamically downscaled GCMs indicate increased extreme precipitation across the western US (Northwest and Southwest NCA regions) for the winter season. Analogously, the dynamically downscaled GCMs both show a general decrease in extreme fall season precipitation across the southeastern US and Mid-Atlantic states. Such consistency, however, is absent for the spring season, wherein the dynamically downscaled GCMs depict changes in projected extreme precipitation of opposite sign across the central US. A one-tailed onesample Wilcoxon test-a nonparametric test that does not require data to be normally distributedis conducted to characterize the robustness of simulated results for each NCA region (table 2). Median projected changes (positive or negative) are statistically significant for every NCA region, except for the dynamically downscaled WRF simulations driven by GFDL RCP 8.5 (2090-2099 with urban expansion over the Southern Great Plains. The signs of significance changes are generally consistent between dynamically downscaled simulations (e.g. across the Midwest, Southwest and Southeast regions). However, for some regions and during some seasons, significant changes of opposite sign are simulated, demonstrating the importance of characterizing impacts from different models. For example, while both dynamically downscaled GCMs reveal statistically significant differences for the Great Plains South region, the sign of the simulated difference is opposing for the winter and spring seasons, with the dynamically downscaled GFDL model simulating a decrease in end-of-century extreme precipitation for these seasons. Conversely, the dynamically downscaled CESM model simulates a statistically significant increase in end-of-century winter and spring extreme precipitation.
Since projected changes in extreme precipitation across NCA regions are generally positive (i.e. extreme precipitation is projected to increase), it is important to quantify the change in the number of days that meet or exceeds the 99th percentile daily total Control precipitation. Figure 5 shows the change in the number of days with extreme precipitation for end relative to start of century for both dynamically downscaled GCMs that include urban expansion. Both dynamically downscaled GCMs indicate an increase in the number of days with extreme precipitation across the western US (Northwest and Southwest NCA regions) for the winter season, which is the largest seasonal contributor to the annual sum. This result is consistent with the increase in extreme precipitation found in these geographical regions during the winter season (e.g. table 2). A north-south gradient is evident, with magnitude increasing northward, ranging from about 10 to 15 additional extreme precipitation winter days across portions of Arizona's Mogollon Rim to more than 30 such additional winter days across portions of the intermountain west (see also table 3). For the winter season, both dynamically downscaled GCMs indicate an increase in the number of extreme days across the western US, a relatively fewer number of additional extreme precipitation days across the eastern US, but no anticipated change across the central US. Although there is broad geographical similarity in the spatial organization of simulated annually averaged extreme precipitation days between the pair of dynamically downscaled GCMs (bottom row of figure 5), the winter and summer season are prominent in highlighting this consistency. Table 3. WRF-simulated change in the number of days-in 10 years (or 10 seasons)-that meet the 99th percentile precipitation exceedance threshold: shown are differences between dynamically downscaled WRF simulations driven by GFDL RCP 8.5 (2090-2099 with urban expansion and Control (2000), and CESM RCP 8.5 (2090-2099) with urban expansion and Control (2000. Units are number of days for the simulated 10 year (or 10 season) period. Asterisk indicates the true possibility of an extreme day for one NCA region at any given day is greater than the reciprocal of corresponding number of days over 10 years (or 10 seasons), based on one-tailed binomial test. All tests are controlled with familywise type-I error rate less than 0.05 (p-value < 0.0007 for each test under the Bonferroni correction for multiple hypothesis tests). Inconsistent statistical significance in the possibility of projected changes in extreme days between CESMRCP85-Urban2100 and GFDLRCP85-Urban2100 are highlighted in bold.

Winter
Spring 2.6 1.9 12.8 * GFDLRCP85-Urban2100 3.4 5.5 * 2.2 3.0 14.1 * Spatial differences between the dynamically downscaled GCMs are also apparent. The dynamically downscaled CESM model indicates a peak in the number of anticipated extreme precipitation days across midwestern states for the spring season while the GFDL model has its peak across the southeast during this time of year. For the fall season, the GFDL model depicts substantial increases in southwestern US extreme precipitation, a feature that is largely absent for the dynamically downscaled CESM model (table 3); however, both dynamically downscaled models simulate little change in extreme precipitation across Southeastern and Mid-Atlantic states. Outside of the winter and, somewhat surprisingly, summer season, compensating seasonal effects are partly responsible for the similarity in simulated annually averaged changes in extreme precipitation between the pair of dynamically downscaled GCMs. Differences in extreme precipitation for dynamically downscaled simulations driven by end-ofcentury GFDL and CESM RCP 8.5, without incorporation of urban expansion, and Control display a similar pattern of extreme precipitation change as the simulations accounting for urban expansion (compare figure 4 to supp. figure 15 (available online at stacks.iop.org/ERL/16/044001/mmedia))). These results imply that projected changes in extreme precipitation are largely driven by emissions of LLGHGs.
We next examine the relationship between urban development, using ∆UHI as a proxy indicator of thermal forcing, and extreme precipitation. We present the relationship between ∆UHI and extreme precipitation change as four quadrant scatter plots that exhibit the following metropolitan region impacts depending on which quadrant an individual metro-region falls into: (a) quadrant I-undergoes an increase in UHI intensity and precipitation; (b) quadrant II-undergoes an increase in UHI intensity and decrease in precipitation; (c) quadrant III-undergoes a decrease in UHI intensity and precipitation; (d) quadrant IV-undergoes a decrease in UHI intensity and increase in precipitation.
Figure 6(a) shows the annual daily averaged change in UHI plotted against the 90th percentile of precipitation. The numbers shown in each of the quadrants indicate the fraction of cities falling in that particular quadrant relative to all cities undergoing change beyond the specified UHI (0.1 K) and precipitation (0.1 mm d −1 ) thresholds. Of those metro-regions undergoing a change greater than our imposed minimum threshold, 50% experience a simulated increase in UHI and 90th percentile of precipitation (quadrant I), while 50% undergo a simulated increase in UHI and a decrease in the 90th percentile of precipitation (quadrant II). We also show the effect of urban development on precipitation for the 95th and 99th percentiles of annually averaged precipitation (figures 6(b) and (c)). For each of the extreme precipitation metrics, the preponderance of cities undergoing changes in extreme precipitation are associated with increased UHI. Metropolitan regions largely fall into quadrant I or II. The imposition of our thresholds do not alter the conclusions drawn, namely that changes in simulated extreme precipitation (either increasing or decreasing precipitation) are broadly associated with an increased UHI intensity. When accounting for increased emissions of LLGHGs we note a shift from primarily quadrants I and II to quadrants I and IV (figures 6(d)-(f)). In agreement with previous single event-scale work our results demonstrate, on a climatic time scale, that while the urban-induced signal can either increase or decrease local to regional-scale precipitation, incorporation of greenhouse gases considerably displaces the probability spectrum toward a general enhancement of extremes. We further note that this shift is consistent between the pair of dynamically downscaled GCM models utilized (table 4).
Although we believe these results are robust at the horizontal resolution utilized, two points of clarification are needed. First, the 20 km grid spacing utilized in our numerical experiments may be insufficient to resolve and ascertain a robust urban-induced drag effect. While the thermal urban signature is well represented at the current resolution (Krayenhoff et al 2018, Broadbent et al 2020, convection-resolving simulations may be necessary to appropriately characterize the local urban-induced trigger for convective systems. In particular, resolving the morphological structure (horizontally and vertically) of urban features that are responsible for UHI-induced convergence zone(s) may be a prerequisite-and requires future attention-for appropriate simulation of convectively driven urban-induced precipitation (e.g. grid spacing of 1-2 km; e.g. Ntelekos et al 2008). Second, prior work has revealed the importance of conditional sampling of a larger subset of characteristic event samples. For example, the New York City UHI has been shown to initiate convective activity over the urban-induced convergence zone resulting from the city's UHI; this effect was apparent during calm regional-scale flows (Bornstein and LeRoy 1990) and has been noted for other cities (Bornstein and Lin 2000). Confirming the climatological fingerprint of urban-induced extreme precipitation may require synoptic typing to retrieve a sample of representative events. Therefore, our results underscore the need for place-based climatological convection-permitting simulations that better account for the urban-associated thermal and drag effects. These physical characteristics are responsible for built environment induced modulation of convergence and divergence patterns, which can lead to initiation of convective systems or, alternatively, may serve as steering mechanisms for already existing convective cells.

Discussion and conclusions
We have conducted a suite of decadal scale WRF simulations examining projected end-of-century annual and seasonal precipitation change across CONUS and the corresponding contribution from the twin forcing agents of urbanization and LLGHGs. Our results indicate increased mean and extreme precipitation for all NCA regions for the concluding relative to the initial decade of the current century. Projected changes in extreme precipitation are greater for the GFDL relative to the CESM dynamically downscaled GCM across all NCA regions except the Southeast. There is broad consistency between the pair of dynamically downscaled GCMs, which show an increase in the number of annual extreme precipitation days across the western US (Northwest and Southwest NCA regions), consistent with the annual increase in extreme precipitation found across these geographical regions. Our results confirm the importance of urban-induced forcing on extreme precipitation for metropolitan areas across the US. We show that while the urban-induced signal can amplify or dampen local to regional-scale precipitation, explicit accounting of LLGHGs displaces the probability spectrum toward a broad enhancement of extremes. The shift toward more extreme precipitation across future CONUS metropolitan regions is consistent between the pair of dynamically downscaled GCM models utilized and highlights the need for improved management policies that address anticipated flooding challenges within urban watersheds. Nevertheless, we affirm the importance of future work undertaking a coordinated and systematically designed suite of place-based, climatological, convection-resolving simulations that better represent the built environment induced dynamical and thermodynamical forcing necessary for appropriate characterization of convective elements. Such coordinated research has the capacity to lead to improved process-based understanding of the key physical mechanisms and feedbacks responsible for extreme precipitation enhancement across urban environments. Importantly, robust results drawn from systematically designed convection-resolving simulations across multiple cities can beneficially guide urban planning across CONUS and global regions where urbanization is and will be occurring at unprecedented scales for decades to come. This is especially important in light of our results, which reveal statistically significant extreme precipitation changes of opposite sign for some NCA regions, highlighting notable simulated differences between the dynamically downscaled GCMs.
In addition, further research with implications for broader urban environmental well-being is needed to facilitate development of resilient cities capable of tackling a diverse array of environmental pressures. The incorporation of biophysical thermal reduction strategies, such as cool and green roofs, can further modify the impacts described herein, and are the subject of ongoing research. Current urban adaptation solutions, as currently designed, are too narrow in their focus, targeting specific outcomes (e.g. reduction in temperature) at the potential expense of others (e.g. air quality), leading to unintended consequences (Fallmann et al 2016, Sailor et al 2016. Optimizing desired outcomes requires a holistic placebased understanding that accounts for large-scale and urban-induced environmental effects, including hydroclimatic (e.g. thermal and precipitation) and air quality impacts contextualized under a socialecological-technological vulnerability framework (Markolf et al 2018). Understanding place-based pressures and vulnerabilities is a necessary condition. Development of wedge-type approaches (Zhao et al 2017) assessing multiple strategies with adaptation alternatives targeting a spectrum of urban-induced but locally contextualized environmental pressures will require participatory engagement, deliberation and assessment of a broader set of future scenarios that highlight tradeoffs among a spectrum of sustainability-oriented indices (Iwaniec et al 2020).

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.