Idealized Simulations of City-Storm Interactions in a Two-Dimensional Framework

Numerous studies have identified spatial variability in convective parameters such as rainfall totals and lightning flashes in the vicinity of large urban areas, yet many questions remain regarding the storm-scale processes that are altered during interaction with a city as well as which urban features are most responsible for storm modification. This study uses an idealized, two-dimensional cloud model to investigate structural and evolutionary changes in a squall line as it passes over a simplified representation of a large city. A parameter space exploration is done in which the parameters of the city—surface temperature and surface roughness length—are systemically increased relative to the region surrounding the idealized city. The resultant suite of simulations demonstrates that storm parameters such as vertical velocity, hydrometeor mass, upward mass flux, and buoyant accelerations are enhanced when the storm passes over the idealized city. No such enhancement occurs in the control simulation without an idealized city.


Introduction
Although comprising only a small percentage of the world's land use, urban areas contain the majority of the world's population. In the United States, 70% of the population currently resides in urbanized areas [1] and it is projected that 68% of the world's population will live in urban areas by 2050 [2]. Convective events that occur over dense, populous urban areas may expose a large number of people to hazards such as flash flooding, damaging wind, large hail, and tornadoes. For this reason, timely and accurate threat information must be communicated to residents of large cities. However, urban environments often present unique challenges to convective forecasters. Large urban areas have been shown to alter the spatial distribution and intensity of convective rainfall around them [3][4][5][6][7][8][9][10][11], with the most commonly observed signal being enhanced rainfall downwind of the urban core. Urban-induced convective initiation is most commonly associated with relatively benign synoptic patterns [12][13][14][15][16] with the resultant thunderstorms typically being short-lived and posing minimal risk for development of convective hazards. However, there is increasing evidence that urban areas can impact organized storm systems and modify the distribution of convective hazards associated with those storms. Huff and Changnon [4] noted that, beginning in the late 1940s, the downwind region of St. Louis, MO experienced more frequent hail events relative to areas upwind and within St. Louis. Reames and Stensrud [17] found that the near-surface rotational and updraft characteristics of supercells are impacted by their proximity to a large urban area. Naylor and Sexton [18] analyzed National Weather Service warning products around several large metropolitan areas within the Midwestern United States. They found that, in most of the cities studied, storm-based severe thunderstorm and tornado warnings were more common on the downwind side of a particular city compared to the upwind side.
There also remains considerable uncertainty regarding which urban features are responsible for storm modification. Huff and Changnon [19] proposed four potential urban processes that may modify convective storms: (1) A thermal mechanism due to the presence of the urban heat island (UHI); (2) an increase in low-level turbulence due to flow obstruction; (3) microphysical modifications due to urban aerosol emissions; and (4) modification of atmospheric moisture content. Various modeling and observational studies have attributed each of these mechanisms to convective variability in some way. Low-level convergence zones induced by the UHI have been associated with areas of frequent convective initiation [20,21]. Surface roughness gradients have been attributed to convergence zones and have been linked to storm bifurcation [22] and downwind enhancement of convective precipitation [23]. Urban aerosols have been shown to play a substantial role in urban storm modification [24][25][26], while spatial gradients in surface moisture over urban areas have also been shown to have an impact [27].
In a real-world setting, it is difficult to isolate these various mechanisms. Idealized storm-scale models allow for controlled variability in the parameters of interest and have been successfully used in the past to investigate the behavior of deep convective storms in particular environments [28][29][30][31][32][33]. The purpose of this study is to use an idealized cloud model to examine the structural changes that occur within an organized convective system as it interacts with a simplified city and to determine if systemically altering characteristics of the city (surface roughness length, urban surface temperature) yields predictable structural changes.

Experiments
Version 18 of the First-Generation Pennsylvania State University-National Center for Atmospheric Research Cloud Model (CM1) [34] is used to conduct a suite of idealized two-dimensional simulations in which a squall line interacts with a simplified representation of an urban area. CM1 is a highly scalable, energy-conserving cloud model specifically designed for idealized simulations of deep convection. The model is configured with a 600 km horizontal domain and a model top of 19.75 km. Stretched grids are used in both the horizontal and vertical directions. The horizontal grid spacing is 250 m within the innermost 120 km of the domain and gradually increases to 1 km outside this region. The vertical grid has 80 total levels with 50 m spacing below 3 km and 500 m spacing above 8.75 km. The model time step is adaptive and automatically adjusts to maintain computational stability. In all simulations the model time step varies between 3-4 s with 8-10 acoustic time steps for each large time step. The lateral boundary conditions are open-radiative, the bottom boundary condition is semi-slip, and the top boundary condition is a rigid wall. A Rayleigh damper is applied above z = 15 km in order to dissipate energy and prevent reflections off the top boundary. The Thompson double-moment microphysics parameterization [35] is used to represent precipitation processes. Subgrid turbulence is represented using the turbulent kinetic energy (TKE) based parameterization and the NASA Goddard parameterization is used for radiative processes.
All simulations utilize a horizontally homogenous base-state environment represented by the Weisman-Klemp analytic thermodynamic profile [28] combined with a unidirectional vertical wind profile from Rotunno et al. [30]. A vertical wind shear of 17.5 m s −1 is confined to the lowest 2.5 km. Convection is initiated at 150 km from the leftmost domain boundary using a thermal bubble. The thermal bubble has a radius of 10 km, is centered at 1.4 km above the surface, and has a maximum potential temperature perturbation of 1 K. The control simulation has a surface skin temperature of 300 K and surface characteristics (albedo, emissivity, and roughness length) of an herbaceous area. The surface roughness length (z 0 ) is set to 0.2 m. In the comparison simulations, a 20-km long region is inserted in the center of the domain to represent the idealized city. The skin temperature inside this region is enhanced to produce a feature comparable to an urban heat island. Surface roughness length is also increased to a value representative of a large urban area. Skin temperature is increased to values of 303, 305, and 307 K-producing differences relative to the area outside the city of 3, 5, and 7 K, respectively. To produce a range of horizontal z 0 gradients, values over the idealized city are increased to 0.5, 1.0, 1.5, and 2 m, respectively. Table 1 shows the values tested within the parameter space. Simulations will be referred to by the magnitude of the difference in skin temperature between the idealized city and outlying area (∆T) as well as the z 0 value of the city in the center of the domain. For example, "5Kz2.0" represents the simulation in which skin temperature within the city is 5 K larger than the surrounding area and the roughness length of the city is set to 2 m. Table 1. Summary of the experiments within the parameter space. All experiments are identical except for two parameters: The surface roughness length (z 0 ) within the idealized city and the skin temperature perturbation (relative to the surrounding area) within the city (∆T).

Trajectory Analysis
In each simulation, 400,000 passive tracers are released into the flow at the beginning of the simulation. These tracers are advected at the model timestep using a 5th order Runge-Kutta advection scheme. The initial location of all tracers is confined between 200 and 600 km in the horizontal domain and between 0.025 and 1.475 km above the surface. Tracers within the storm updraft at a specific analysis time are followed backwards in time to determine their previous positions with a temporal resolution of 1 min. Model-calculated fields (e.g., buoyancy, vertical velocity) are interpolated to the tracer positions, also with a temporal resolution of 1 min.

Results
All simulations produce a squall line that develops from the initial thermal perturbation and travels left-to-right across the domain. By t = 120 min, the squall line has not yet passed over the idealized city in any of the simulations, but differences in the low-level thermodynamic structure over the idealized cities are readily apparent (Figure 1). In the control simulation (0Kz0.2), temperature is horizontally homogeneous in the environment ahead of the squall line. The same is true in all ∆T = 0 K simulations. However, as ∆T increases, a dome-like structure in the near-surface temperature develops (see 3Kz0.2, 5Kz0.2, 7Kz0.2 simulations in Figure 1). As expected, the magnitude of the temperature excess in this thermal dome increases as ∆T is increased. Additionally, in simulations with ∆T > 0 K, increasing z 0 results in a stronger and deeper thermal dome. For example, in the 7Kz0.2 simulation, the thermal dome structure does not extend above z = 0.5 km and temperature values within the dome are approximately 1 K warmer than areas outside the dome. In the 7Kz2.0 simulation, the thermal perturbation at the center of the domain extends up to 0.75 km and temperature values are 3-4 K larger than outside the dome. The depth of this dome is confined below z = 1 km in all simulations, in agreement with observations of real urban heat islands [36]. A further investigation (not shown) revealed that the simulations with enhanced z 0 produced larger eddy mixing coefficients, which more efficiently mixed the warmer temperatures from the heated surface. Figure 2 shows the vertical structure of potential temperature (θ) at t = 120 min. In the control simulation (0Kz0.2), a low-level stratification is evident with θ decreasing with height across the domain. This same structure is also present in the other 0 K simulations as z 0 is increased over the city (0Kz0.5, 0Kz1.0, 0Kz1.5, 0Kz2.0). The introduction of a heated surface disrupts this stratification directly over the idealized city, leading to instances of constant θ with height in the lowest kilometer above the surface (3Kz0.2, 5Kz0.2, 7Kz0.2). This feature agrees with observations of real urban heat islands [37] although the depth of the neutral stratification varies across the parameter space simulations. As previously discussed, vertical mixing is enhanced when z 0 is increased. In the simulations with enhanced ∆T, Atmosphere 2020, 11, 707 4 of 24 increasing the z 0 gradient results in a deeper layer with near neutral θ stratification. For a given ∆T, the deepest heat islands are found with the largest z 0 values.
Atmosphere 2020, 10, x FOR PEER REVIEW 4 of 24 Figure 1. Cross-section of temperature (K) at t = 120 min. Each row represents simulations with a particular skin temperature perturbation within the idealized city while each column represents simulations with a particular z0 value within the city. The alterations to the low-level thermodynamic structure within the city have a dramatic impact on instability. Table 2 shows values of surface-based convective available potential energy (CAPE) and convective inhibition (CIN) averaged across the idealized city at t = 120 min. In the control environment, CAPE is 1932 J kg −1 and CIN is −46 J kg −1 . Increasing z0 without increasing ΔT (top row   The alterations to the low-level thermodynamic structure within the city have a dramatic impact on instability. Table 2 shows values of surface-based convective available potential energy (CAPE) and convective inhibition (CIN) averaged across the idealized city at t = 120 min. In the control environment, CAPE is 1932 J kg −1 and CIN is −46 J kg −1 . Increasing z0 without increasing ΔT (top row  The alterations to the low-level thermodynamic structure within the city have a dramatic impact on instability. Table 2 shows values of surface-based convective available potential energy (CAPE) and convective inhibition (CIN) averaged across the idealized city at t = 120 min. In the control environment, CAPE is 1932 J kg −1 and CIN is −46 J kg −1 . Increasing z 0 without increasing ∆T (top row in Table 2) results in a small CAPE increase and a negligible change in CIN. A much more noticeable change occurs when ∆T > 0 K. For a given value of z 0 , increasing ∆T results in a substantial increase in CAPE and a reduction in CIN (see the 3Kz0.2, 5Kz0.5, and 7Kz0.2 simulations). When ∆T > 0 K, increasing the z 0 gradient also results in increased CAPE and decreased CIN (see the 7Kz0.5, 7Kz1.0, 7Kz1.5, and 7Kz2.0 simulations). When comparing the control (0Kz0.2) and 7Kz2.0 simulations, CAPE increases from 1932 to 2855 J kg −1 (47% increase) while CIN decreases from −46 to −1 J kg −1 . Table 2. Surface-based CAPE (J kg −1 ) and CIN (J kg −1 , italicized) over the idealized city at t = 120 min. Variations in CAPE and CIN between simulations should be associated with changes to updraft intensity as storms approach and pass over the city. While this is true in the simulations, the relationship is not monotonic. In other words, the simulations with largest CAPE/lowest CIN do not necessarily yield the strongest updrafts. Table 3 shows the maximum instantaneous updraft velocity during a 60 min window (t = 140 min to t = 200 min) focused on the time period when the squall line passes over the idealized city. In the control simulation (0Kz0.2), the domain maximum vertical velocity is 41 m s −1 -the smallest value within the parameter space. When ∆T = 0 K, inducing a horizontal z 0 gradient results in larger values of maximum vertical velocity. Increasing ∆T while the horizontal z 0 gradient is zero or relatively modest (z 0 = 0.2 m, 0.5 m) also results in larger vertical velocities. The largest vertical velocities are found in the 5Kz0.5 and 7Kz0.5 simulations. Conversely, increasing ∆T in the simulations with larger z 0 gradients (z 0 = 1.5 m, 2.0 m) typically results in a much smaller increase relative to the control. Aside from the control simulation, the smallest instantaneous vertical velocities are found with large z 0 and large ∆T (e.g., the 5Kz2, 7Kz1.0, 7Kz1.5, and 7Kz2.0 simulations). While these simulations have the largest CAPE over the city, they also have the least amount of CIN. In these simulations, convection is initiated downwind of the city and interacts with the main squall line. This interaction may have negative consequences to storm intensity.   Atmosphere 2020, 11, 707 6 of 24 Figure 3 shows the time-series of domain maximum vertical velocity over the duration of each simulation. During the first 80 min, all simulations produce an identical evolution with smaller differences emerging after this time. Each simulation exhibits a sharp ramp-up period over the first 50 min of simulation followed by a quasi-steady period that lasts until approximately t = 120 min. From t = 120 min through t = 180 min, most simulations exhibit a steady increase in maximum vertical velocity to a new quasi-steady state of roughly 35-40 m s −1 . After t = 200 min, most of the simulations show a gradual weakening in the maximum updraft. Despite the overall similarities, differences in evolution are evident. For example, the 0Kz0.2 simulation exhibits updraft weakening between t = 180 min and t = 210 min followed by a restrengthening. In simulations with ∆T = 0 K and enhanced z 0 , the weakening is not present and the maximum vertical velocity is enhanced. For a given value of z 0 , an increase in ∆T does not appear to have a notable impact on the domain maximum vertical velocity in most simulations. The exception to this trend is the z 0 = 0.5 m experiments. The 5Kz0.5 and 7Kz0.5 cases have an abrupt spike in maximum vertical velocity centered on the time of interaction with the idealized city that is not present in the 0Kz0.5 simulation. The 3Kz0.5 simulation does exhibit a spike in vertical velocity, however it occurs approximately 20-30 min after interaction with the city. energy of upward motion within the domain. A term herein referred to as updraft kinetic energy (UKE) is defined as: where mi is the mass within the ith grid cell and wi is the vertical velocity within the ith cell. The summation is over all points where w > 1 m s −1 and within the cloud (total hydrometeor mixing ratio > 1 g kg −1 ). The time-series of UKE for all simulations are shown in Figure 4. In the control simulation (0Kz0.2), UKE fluctuates below 2 × 10 6 J until approximately t = 200 min when it increases to a maximum of over 4 × 10 6 J between 250-300 min. When the z0 gradient is increased and ΔT = 0 K, the UKE maximum occurs earlier in time. For example, in the 0Kz2 simulation, UKE exceeds 4 × 10 6 J before t = 200 min, approximately 60 min earlier than in the control simulation. When ΔT is increased but z0 remains the same as the surroundings (e.g., 3Kz0.2, 5Kz0.2, 7Kz0.2), a similar trend is seen, with the maximum at the end of the simulation moving to progressively earlier times. Other differences in the UKE evolution are also evident, but they do not exhibit a consistent trend. In the 3Kz0.2 simulation, there is a sharp increase in UKE to over 3 × 10 6 J during interaction with the idealized city. When ΔT is increased further (5Kz0.2, 7Kz0.2), that increase is no longer present. In simulations where both z0 and ΔT are increased, a substantial increase in UKE often (but not always) occurs during or shortly after the city-storm interaction period (see the 3Kz0.5, 3Kz2.0, 5Kz0.5, 5Kz1.0, 5Kz1.5, 7Kz0.5, and 7Kz1.0 simulations). The largest increase is seen in the 5Kz0.5 simulation. In other simulations (3Kz1.5, 7Kz1.5, 7Kz2.0), no such increase is apparent.  More distinct changes in the updraft structure are evident when evaluating the total kinetic energy of upward motion within the domain. A term herein referred to as updraft kinetic energy (UKE) is defined as: where m i is the mass within the ith grid cell and w i is the vertical velocity within the ith cell. The summation is over all points where w > 1 m s −1 and within the cloud (total hydrometeor mixing ratio > 1 g kg −1 ). The time-series of UKE for all simulations are shown in Figure 4. In the control simulation (0Kz0.2), UKE fluctuates below 2 × 10 6 J until approximately t = 200 min when it increases to a maximum of over 4 × 10 6 J between 250-300 min. When the z 0 gradient is increased and ∆T = 0 K, the UKE maximum occurs earlier in time. For example, in the 0Kz2 simulation, UKE exceeds 4 × 10 6 J before t = 200 min, approximately 60 min earlier than in the control simulation. When ∆T is increased but z 0 remains the same as the surroundings (e.g., 3Kz0.2, 5Kz0.2, 7Kz0.2), a similar trend is seen, with Atmosphere 2020, 11, 707 7 of 24 the maximum at the end of the simulation moving to progressively earlier times. Other differences in the UKE evolution are also evident, but they do not exhibit a consistent trend. In the 3Kz0.2 simulation, there is a sharp increase in UKE to over 3 × 10 6 J during interaction with the idealized city. When ∆T is increased further (5Kz0.2, 7Kz0.2), that increase is no longer present. In simulations where both z 0 and ∆T are increased, a substantial increase in UKE often (but not always) occurs during or shortly after the city-storm interaction period (see the 3Kz0.5, 3Kz2.0, 5Kz0.5, 5Kz1.0, 5Kz1.5, 7Kz0.5, and 7Kz1.0 simulations). The largest increase is seen in the 5Kz0.5 simulation. In other simulations (3Kz1.5, 7Kz1.5, 7Kz2.0), no such increase is apparent.   Figure 6 shows the vertical velocity structure 40 min later, after the storms have passed over the idealized city. In most simulations, the primary updraft is stronger and extends to a higher altitude relative to the structure at t = 130 min. In addition to stronger updrafts, the cloud area has also expanded in all cases. In the cases with large ∆T and large z 0 gradients (e.g., 5Kz2.0, 7Kz1.0, 7Kz1.5, 7Kz2.0), the convection on the downwind side of the city has become more expansive compared to the earlier time, with cells extending up to 70 km downwind of the city center. In the absence of a heat island (∆T = 0 K), increasing the horizontal z 0 gradient results in stronger updrafts. Compared to the control simulation, the 0Kz0.5, 0Kz1.0, and 0Kz2.0 have stronger updrafts at t = 170 min. In the z 0 = 0.2 m simulations, increasing ∆T also results in stronger and deeper updrafts. When compared to the control simulation, the 3Kz0.2, 5Kz0.2, and 7Kz0.2 simulations all have stronger updrafts, particularly between 5 and 10 km above the surface. When both ∆T and the z 0 gradient are increased, most simulations contain updrafts that are stronger and deeper compared to the control simulation, however systematically increasing one parameter (either ∆T or z 0 ) does not always lead to stronger updrafts. For example, the updraft in the 3Kz0.2 simulation appears stronger than in the 3Kz0.5 or 3K2.0 simulations. The strongest updrafts at t = 170 min occur in the 0Kz2.0, 5Kz0.5, 5Kz1.0, and 7Kz1.0 simulations.  Changes in the vertical velocity before and after interaction with the idealized city are further emphasized in Figure 7. Changes to the updraft intensity and/or size are smallest in magnitude in the control simulation (0Kz0.2). As the surface z0 gradient is increased (top row of Figure 7), a more noticeable updraft strengthening is apparent. When ΔT is increased but the z0 gradient remains zero   Changes in the vertical velocity before and after interaction with the idealized city are further emphasized in Figure 7. Changes to the updraft intensity and/or size are smallest in magnitude in the control simulation (0Kz0.2). As the surface z0 gradient is increased (top row of Figure 7), a more noticeable updraft strengthening is apparent. When ΔT is increased but the z0 gradient remains zero (0Kz0.2, 3Kz0.2, 5Kz0.2, and 7Kz0.2 simulations), updraft strengthening is evident at mid to upper levels. When both ΔT and the surface z0 gradient are enhanced, changes to the updraft structure range from minimal (e.g., 3Kz0.5) to expansive (e.g., 5Kz0.5). Changes in the vertical velocity before and after interaction with the idealized city are further emphasized in Figure 7. Changes to the updraft intensity and/or size are smallest in magnitude in the control simulation (0Kz0.2). As the surface z 0 gradient is increased (top row of Figure 7), a more noticeable updraft strengthening is apparent. When ∆T is increased but the z 0 gradient remains zero (0Kz0.2, 3Kz0.2, 5Kz0.2, and 7Kz0.2 simulations), updraft strengthening is evident at mid to upper levels. When both ∆T and the surface z 0 gradient are enhanced, changes to the updraft structure range from minimal (e.g., 3Kz0.5) to expansive (e.g., 5Kz0.5). The blue contour represents the hydrometeor mixing ratio of 0.02 g kg −1 at t=170 min. The colored field was calculated by subtracting the vertical velocity field at t = 170 min from the same field at t = 130 min after adjusting for a storm motion in the earlier field using the mean environmental winds.
Variability in vertical velocity between simulations is associated with changes to low-level divergence. Figure 8 shows divergence and vertical velocity at t = 145 min in all simulations. The impact of horizontal z0 gradients is most evident in the ΔT = 0 K simulations (top row of Figure 8). For example, in the 0Kz0.2 simulation there is a small area of divergence centered on the location of the idealized city (x = 0 km). As the z0 gradient is increased, convergence increases in both magnitude The blue contour represents the hydrometeor mixing ratio of 0.02 g kg −1 at t=170 min. The colored field was calculated by subtracting the vertical velocity field at t = 170 min from the same field at t = 130 min after adjusting for a storm motion in the earlier field using the mean environmental winds.
Variability in vertical velocity between simulations is associated with changes to low-level divergence. Figure 8 shows divergence and vertical velocity at t = 145 min in all simulations. The impact of horizontal z 0 gradients is most evident in the ∆T = 0 K simulations (top row of Figure 8). For example, in the 0Kz0.2 simulation there is a small area of divergence centered on the location of the idealized city (x = 0 km). As the z 0 gradient is increased, convergence increases in both magnitude and area over the city. In the 0Kz2.0 simulation, convergence extends almost to the lateral edges of the city (x = −10 km to x = 10 km) and extends upward beyond 2 km above the surface. This area of convergence is associated with updrafts directly over the city above z = 1 km. While the simulations with ∆T > 0 K have somewhat chaotic patterns in convergence/divergence (rows 2-4 in Figure 8), there are coherent patterns in relation to vertical velocity. Increasing ∆T results in areas of updraft extending closer to the surface (see the 0Kz0.2, 3Kz0.2, and 5Kz0.2 simulations). In general, low-level divergence is associated with downdrafts while areas of low-level convergence are associated with updrafts.
Atmosphere 2020, 10, x FOR PEER REVIEW 10 of 24 Low-level divergence patterns in three of the ΔT = 0 K simulations (0Kz0.2, 0Kz1.0, and 0Kz2.0) at t = 145 min are shown in Figure 9. At this time, the leading edge of the squall line is interacting with the leftmost edge of the idealized city (x = -10 km). As the storm approaches, outflow from the negative-x direction is directed over the city. On the opposite side of the city, low-level inflow approaches the leading edge of the storm from the positive-x direction. At the interface of the inflow and outflow, convergence is present. In the 0Kz0.2 simulation, the near-surface convergence is fairly weak and there is a small area of divergence directly over the idealized city (x = 0 km). When the surface z0 gradient is increased, the inflow into the storm is slowed over the idealized city, resulting in a broader and stronger area of near-surface convergence in the 0Kz1.0 and 0Kz2.0 simulations.
Updraft strengthening is also associated with changes in buoyant accelerations. Buoyancy (B) is calculated as , where g is gravitational acceleration, is density potential temperature, and subscript 0 indicates the base-state value. Buoyancy values along forward-integrated parcel trajectories are analyzed and traced backwards for 50 min. Trajectories are filtered based on a vertical velocity threshold at the defined analysis time. For example, Figure 10 only shows the location history of parcels with a vertical velocity greater than 1 m s −1 at t = 120 min. These parcels are traced backwards for 50 min to determine their origin location and buoyancy history. In all simulations, the majority of parcels originate ahead of the storm, although some parcels do appear to originate in the storm outflow. In all simulations, the updraft parcels exhibit a midlevel maximum in B between 4 Low-level divergence patterns in three of the ∆T = 0 K simulations (0Kz0.2, 0Kz1.0, and 0Kz2.0) at t = 145 min are shown in Figure 9. At this time, the leading edge of the squall line is interacting with the leftmost edge of the idealized city (x = -10 km). As the storm approaches, outflow from the negative-x direction is directed over the city. On the opposite side of the city, low-level inflow approaches the leading edge of the storm from the positive-x direction. At the interface of the inflow and outflow, convergence is present. In the 0Kz0.2 simulation, the near-surface convergence is fairly weak and there is a small area of divergence directly over the idealized city (x = 0 km). When the surface z 0 gradient is increased, the inflow into the storm is slowed over the idealized city, resulting in a broader and stronger area of near-surface convergence in the 0Kz1.0 and 0Kz2.0 simulations.
Updraft strengthening is also associated with changes in buoyant accelerations. Buoyancy (B) is calculated as = g θ ρ −θ ρ0 θ ρ0 , where g is gravitational acceleration, θ ρ is density potential temperature, and subscript 0 indicates the base-state value. Buoyancy values along forward-integrated parcel trajectories are analyzed and traced backwards for 50 min. Trajectories are filtered based on a vertical velocity threshold at the defined analysis time. For example, Figure 10 only shows the location history of parcels with a vertical velocity greater than 1 m s −1 at t = 120 min. These parcels are traced backwards for 50 min to determine their origin location and buoyancy history. In all simulations, the majority of parcels originate ahead of the storm, although some parcels do appear to originate in the storm outflow. In all simulations, the updraft parcels exhibit a midlevel maximum in B between 4 and 10 km above the surface. The magnitude of the maximum B values is similar between simulations. The same is true for the spatial extent of the maximum. One of the more notable differences is seen near the surface. As ∆T and the z 0 gradient are increased, near surface B values also increase. By t = 165 min (after passing over the idealized city) there are much more apparent differences between simulations (Figure 11). In all simulations, the overall area covered by the largest B values has increased both horizontally and vertically. For example, at t = 120 min (Figure 10) in the 3Kz0.2 simulation, B values greater than 1 × 10 −4 m s −2 are concentrated mainly in a narrow horizontal area between 5 and 7.5 km above the surface. By t = 165 min (Figure 11), B values above this threshold are By t = 165 min (after passing over the idealized city) there are much more apparent differences between simulations (Figure 11). In all simulations, the overall area covered by the largest B values has increased both horizontally and vertically. For example, at t = 120 min ( Figure 10) in the 3Kz0.2 simulation, B values greater than 1 × 10 −4 m s −2 are concentrated mainly in a narrow horizontal area between 5 and 7.5 km above the surface. By t = 165 min (Figure 11), B values above this threshold are found below 4 and above 10 km. The horizontal extent also increased substantially. While all simulations show an increase in B between 120 and 165 min, the increase is greatest in the simulations with larger ∆T and larger z 0 gradients. The smallest differences in B are found in the 0Kz0.2, 0Kz0.5, 3Kz0.2, and 3Kz0.5 simulations. Updrafts within the cells initiated downwind of the idealized city experience substantial buoyant accelerations in the 7Kz1.0, 7Kz1.5, and 7Kz2.0 simulations.
Atmosphere 2020, 10, x FOR PEER REVIEW 12 of 24 Figure 10. Time-space evolution of forward-integrated parcel trajectories that are located in the updraft at t = 120 min. The position of each parcel over the previous 50 min is shown, with a temporal resolution of 1 min. At each location, the parcel's location marker is colored by B. The black line is the 0.02 g kg −1 contour hydrometeor mixing ratio at t = 120 min. Figure 10. Time-space evolution of forward-integrated parcel trajectories that are located in the updraft at t = 120 min. The position of each parcel over the previous 50 min is shown, with a temporal resolution of 1 min. At each location, the parcel's location marker is colored by B. The black line is the 0.02 g kg −1 contour hydrometeor mixing ratio at t = 120 min. Figure 10. Time-space evolution of forward-integrated parcel trajectories that are located in the updraft at t = 120 min. The position of each parcel over the previous 50 min is shown, with a temporal resolution of 1 min. At each location, the parcel's location marker is colored by B. The black line is the 0.02 g kg −1 contour hydrometeor mixing ratio at t = 120 min.  Further insights into the temporal evolution of updraft accelerations can be gained by analyzing contoured frequency by altitude diagrams (CFADs; Yuter and Houze [38]) of B. Figure 12 shows CFADs of maximum B along parcel trajectories for all simulations at t = 120 min. At this time, CFADs of B are very similar. In all simulations, most parcels achieve maximum buoyancy below 8 km and all parcels achieve a maximum B less than 2.5 × 10 −4 m s −2 .
achieve maximum B above 8 km, but all maximum values of B in the control simulation remain less than 2.5 × 10 −4 m s −2 . In contrast, updrafts in simulations with larger ΔT and/or larger z0 gradient vary from the control in several important ways. First, parcel trajectories that achieve a large positive B (greater than 2.0 × 10 −4 m s −2 ) become more common compared to the earlier time and parcels achieve this large B over a more diverse range of heights. Second, the maximum B of any one parcel increases. In many of the simulations, some parcels achieve a maximum B greater than 3 × 10 −4 m s −2 . Figure 12. CFADs of maximum buoyancy along forward-integrated trajectory parcels located within an updraft at t = 120 min for each simulation. Maximum buoyancy is determined by examining the parcel's history over the previous 50 min. Only parcels that achieve positive buoyancy are shown. Figure 12. CFADs of maximum buoyancy along forward-integrated trajectory parcels located within an updraft at t = 120 min for each simulation. Maximum buoyancy is determined by examining the parcel's history over the previous 50 min. Only parcels that achieve positive buoyancy are shown. Figure 13 shows CFADs of B at t = 165 min. In the control simulation (0Kz0.2), maximum B values are similar to those at t = 120 min ( Figure 12). One difference at t = 165 min is that some parcels achieve maximum B above 8 km, but all maximum values of B in the control simulation remain less than 2.5 × 10 −4 m s −2 . In contrast, updrafts in simulations with larger ∆T and/or larger z 0 gradient vary from the control in several important ways. First, parcel trajectories that achieve a large positive B (greater than 2.0 × 10 −4 m s −2 ) become more common compared to the earlier time and parcels achieve this large B over a more diverse range of heights. Second, the maximum B of any one parcel increases. In many of the simulations, some parcels achieve a maximum B greater than 3 × 10 −4 m s −2 .
The vertical perturbation pressure gradient force (vppgf) calculated along each parcel trajectory is shown in Figure 14. The parcels shown in Figure 14 are located in the updraft at t = 165 min and traced backwards for 50 min. All simulations contain a 'patch' of large values of vppgf at low levels. In most of the simulations, a narrow band of larger vppgf extends upwards in the main updraft. One of the larger differences between the simulations is seen near the cloud top, which is above 10 km. It can be seen that for a given ∆T, increasing the horizontal z 0 gradient typically results in a stronger vppgf near the cloud top. The increased vppgf near the cloud top indicates that these simulations contain wider updrafts, as shown in previous studies [39][40][41]. Atmosphere 2020, 10, x FOR PEER REVIEW 14 of 24 Figure 13. Same as Figure 12, except at t = 165 min.
The vertical perturbation pressure gradient force (vppgf) calculated along each parcel trajectory is shown in Figure 14. The parcels shown in Figure 14 are located in the updraft at t = 165 min and traced backwards for 50 min. All simulations contain a 'patch' of large values of vppgf at low levels. In most of the simulations, a narrow band of larger vppgf extends upwards in the main updraft. One of the larger differences between the simulations is seen near the cloud top, which is above 10 km. It can be seen that for a given ΔT, increasing the horizontal z0 gradient typically results in a stronger vppgf near the cloud top. The increased vppgf near the cloud top indicates that these simulations contain wider updrafts, as shown in previous studies [39][40][41].    Figure 11, except the parcel location is shaded by the vertical perturbation pressure gradient force. Figure 15 displays an updraft mass flux (F m ), which is calculated as F m = max w i,k , 0 dx 2 , where w i,k is the vertical velocity at the ith horizontal and kth vertical grid point, and dx is the horizontal grid spacing. At t = 120 min in the control simulation (0Kz0.2), F m increases steadily from the surface until reaching a maximum value around 4 km above the surface at which point F m decreases steadily until approximately 12 km above the surface. Above 12 km, F m exhibits little change with height. At t = 165 min, the F m profile in the control simulation is nearly identical to that at the earlier time. At both t = 120 min and t = 165 min, maximum F m in the control simulation is approximately 2.2 × 10 7 kg s −1 . The primary difference between the two times is a slight increase in F m aloft beginning at approximately 7 km above the surface. In contrast, when ∆T = 0 K, increasing the horizontal z 0 gradient results in progressively larger F m after interaction with the idealized city. The F m maximum also tends to occur at a higher altitude as z 0 is increased. In the 0Kz2.0 simulation, the maximum F m at t = 165 min increases to 3.0 × 10 7 kg s −1 , which is an increase of over 36% compared to the control simulation. Increasing ∆T while holding the horizontal z 0 gradient to zero results in a similar trend. In the 3Kz0.2 simulation, the peak F m value reaches 2.6 × 10 7 kg s −1 between 4 and 6 km above the surface while the 5Kz0.2 simulation has a maximum F m of 2.9 × 10 7 kg s −1 . Increasing ∆T further (7Kz0.2) results in a slightly smaller maximum F m (2.7 × 10 7 kg s −1 ). For a fixed value of ∆T (z 0 ), increasing z 0 (∆T) typically results in increased F m . The largest F m value is seen in the 7Kz2.0 simulation, where the maximum of 3.7 × 10 7 kg s −1 is 68% larger compared to the maximum in the control simulation. Stronger, more vigorous updrafts may increase the likelihood of surface or near-surface convective hazards such as flash flooding, hail, or damaging winds. Figure 16 shows time series plots of domain total graupel/hail mass (mgraup) in each simulation. In the control simulation, there is a steady upward trend from t = 100 min through approximately t = 170 min, with a maximum value of 4.7 × 10 7 kg. After the maximum is achieved, there is a slight reduction in mgraup before the values become steady. Within the tested parameter space, the control contains the smallest maximum mgraup value during the 100 to 200 min time window. When the horizontal z0 gradient is increased but ΔT remains zero, a substantial increase in mgraup (relative to the control simulation) occurs. In the 0Kz0.5 simulation, mgraup surpasses the maximum in the control simulation before the squall line passes over the city. By t = 170 min, after the storm has passed over the city, mgraup exceeds 8 × 10 7 kg. In the 0Kz2.0 simulation, maximum mgraup exceeds 10 × 10 7 kg at t = 200 min. Increasing ΔT while z0 = 0.2 m results in similar increases to mgraup. In the 3Kz0.2 simulation, the time evolution of mgraup is similar to the Stronger, more vigorous updrafts may increase the likelihood of surface or near-surface convective hazards such as flash flooding, hail, or damaging winds. Figure 16 shows time series plots of domain total graupel/hail mass (m graup ) in each simulation. In the control simulation, there is a steady upward trend from t = 100 min through approximately t = 170 min, with a maximum value of 4.7 × 10 7 kg. After the maximum is achieved, there is a slight reduction in m graup before the values become steady. Within the tested parameter space, the control contains the smallest maximum m graup value during the 100 to 200 min time window. When the horizontal z 0 gradient is increased but ∆T remains zero, a substantial increase in m graup (relative to the control simulation) occurs. In the 0Kz0.5 simulation, m graup surpasses the maximum in the control simulation before the squall line passes over the city. By t = 170 min, after the storm has passed over the city, m graup exceeds 8 × 10 7 kg. In the 0Kz2.0 simulation, maximum m graup exceeds 10 × 10 7 kg at t = 200 min. Increasing ∆T while z 0 = 0.2 m results in similar increases to m graup . In the 3Kz0.2 simulation, the time evolution of m graup is similar to the control simulation except that the peak is greater-5.8 × 10 7 kg compared to 4.7 × 10 7 kg. In addition, during interaction with the idealized city the rate of increase (inferred from the slope of the time-series line) is substantially greater in the 3Kz0.2 simulation compared to the control. In the 7Kz0.2 simulation, the rate of increase is similar to the 3Kz0.2 simulation but the period of increase is longer and the peak in m graup is also substantially larger. When both ∆T and the horizontal z 0 gradient are increased, many of the simulations exhibit a drastic increase in m graup during (or slightly after) interaction with the city (see 3Kz1.5, 5Kz0.5, 5Kz1.0, 5Kz1.5, 7Kz0.5, 7Kz1.0, 7Kz1.5). An additional difference can be seen prior to interaction with the idealized city. In the control simulation, m graup remains between 2 × 10 7 kg and 3 × 10 7 kg from t = 100 min through t = 130 min before steadily increasing. Many of the simulations with enhanced ∆T and/or z 0 exhibit a decrease in m graup around this time, with several of the simulations exhibiting rather sharp decreases just prior to interacting with the idealized city (0Kz0.5, 3Kz0.2, 3Kz1.0, 5Kz0.2). A similar evolution is apparent in the domain total rain mass (mrain). Figure 17 shows that in the control simulation, the mrain is relatively constant from t = 100 min through t = 200 min, except for a slight increase in mrain after t = 180 min which quickly subsides. In contrast, when ΔT remains zero, simulations with an enhanced z0 gradient produce more substantial increases in mrain following interaction with the city. Increasing ΔT also tends to result in larger mrain. All of the ΔT = 7 K simulations have a maximum mrain of greater than 4.0 × 10 7 kg during the time window shown in Figure 17, whereas for the ΔT = 0 K and ΔT = 3 K simulations, only 0Kz1.0 and 3Kz1.0 contain mrain exceeding 4.0 × 10 7 kg. In addition, many of the simulations (0Kz2.0, 3Kz1.0, 3Kz1.5, 5Kz0.2, 5Kz0.5, 5Kz1.0, 5K2.0, 7Kz0.2, 7Kz1.0, 7K1.5, 7K2.0) exhibit a pronounced increase in mrain when the squall line is passing over the idealized city.
While nearly every simulation with increased ΔT and/or z0 produces larger mrain relative to the control simulation, there is no predictable pattern in mrain when either ΔT or z0 is increased. For A similar evolution is apparent in the domain total rain mass (m rain ). Figure 17 shows that in the control simulation, the m rain is relatively constant from t = 100 min through t = 200 min, except for a slight increase in m rain after t = 180 min which quickly subsides. In contrast, when ∆T remains zero, simulations with an enhanced z 0 gradient produce more substantial increases in m rain following interaction with the city. Increasing ∆T also tends to result in larger m rain . All of the ∆T = 7 K simulations have a maximum m rain of greater than 4.0 × 10 7 kg during the time window shown in Figure 17, whereas for the ∆T = 0 K and ∆T = 3 K simulations, only 0Kz1.0 and 3Kz1.0 contain m rain exceeding 4.0 × 10 7 kg. In addition, many of the simulations (0Kz2.0, 3Kz1.0, 3Kz1.5, 5Kz0.2, 5Kz0.5, 5Kz1.0, 5K2.0, 7Kz0.2, 7Kz1.0, 7K1.5, 7K2.0) exhibit a pronounced increase in m rain when the squall line is passing over the idealized city.
Atmosphere 2020, 10, x FOR PEER REVIEW 18 of 24 Figure 17. Same as Figure 16, except for rain mass. Figure 18 shows cross-sections of the graupel/hail mixing ratio (qg) after interaction with the idealized city. In the control simulation, the largest qg values are found between 3 and 5 km above the surface, with a maximum value of 6.6 g kg −1 . Increasing z0 over the city while ΔT = 0 K results in larger values of qg as well as an increase in the area of qg > 1 g kg −1 . Increased qg is also found when ΔT is increased while z0 remains 0.2 m. The largest values of qg are found in the simulations where both ΔT and the horizontal z0 gradient are increased. In particular, a maximum qg in both the 5Kz1.0 and 7Kz0.5 simulations exceed 13 g kg −1 . Within the tested parameter space, only the 3Kz0.2 simulation yields a smaller maximum qg value at t = 190 min compared to the control simulation.
Interaction with the idealized city also has an impact on the near-surface winds. Figure 19 shows the time-series of the maximum 10-m wind speed in each simulation. In the control simulation, the maximum fluctuates around 10 m s −1 between t = 100 min and t = 170 min. After 170 min, there is a sharp, brief increase to a peak value of 14 m s −1 followed by a slight decrease. This value (14 m s −1 ) represents the smallest peak 10-m wind speed in any of the simulations. In the 0Kz0.5 simulation, the peak maximum is larger compared to the control simulation (24.5 m s −1 vs. 14 m s −1 ) and occurs approximately 10 min earlier. While this peak value is short-lived, values of roughly 15 m s −1 persist until the end of the simulation (not shown). The remaining ΔT = 0 K simulations also exhibit steady maximum winds of roughly 15 m s −1 after interaction with the city, although none have the drastic peak seen in the 0Kz0.5 simulation. Increasing ΔT also results in increased maximum values. For example, in the 3Kz0.2 simulation, the maximum wind speed of 16 m s −1 occurs 36 min prior to the maximum in the control simulation, when the squall line begins interacting with the idealized city. Further increasing ΔT to 5 K results in even stronger 10-m winds. A value of 16 m s −1 occurs during the interaction with the idealized city, and another peak of 21 m s −1 occurs around t = 200 min. In the 7Kz0.2 simulation, the maximum 10-m wind is larger than in the control and 3Kz0.2, but weaker than in 5Kz0.2. Many of the simulations with enhanced ΔT and horizontal z0 gradients produce sharp peaks in the maximum 10-m wind during or shortly after interaction with the idealized city. The strongest winds are found in the 3Kz1.0 (23.6 m s −1 at 222 min), 3Kz2.0 (25 m s −1 at 203 min), 5Kz0.5 While nearly every simulation with increased ∆T and/or z 0 produces larger m rain relative to the control simulation, there is no predictable pattern in m rain when either ∆T or z 0 is increased. For example, the 3Kz0.2 simulation does not show much of an increase in m rain (relative to the control simulation) during or after interaction with the idealized city, whereas the 3Kz1.0 simulation contains a substantial increase in m rain just prior to t = 200 min. However, further increasing the z 0 gradient (3Kz1.5) results in a smaller peak in m rain , while the 3Kz2.0 does not contain a sharp peak at all. Increasing ∆T also yields inconsistent results. The 5Kz0.2 and 7Kz0.2 have similar increases in m rain between t = 180 min and t = 200 min. The same is true for the 5Kz0.5 and 7Kz0.5 simulations. Figure 18 shows cross-sections of the graupel/hail mixing ratio (q g ) after interaction with the idealized city. In the control simulation, the largest q g values are found between 3 and 5 km above the surface, with a maximum value of 6.6 g kg −1 . Increasing z 0 over the city while ∆T = 0 K results in larger values of q g as well as an increase in the area of q g > 1 g kg −1 . Increased q g is also found when ∆T is increased while z 0 remains 0.2 m. The largest values of q g are found in the simulations where both ∆T and the horizontal z 0 gradient are increased. In particular, a maximum q g in both the 5Kz1.0 and 7Kz0.5 simulations exceed 13 g kg −1 . Within the tested parameter space, only the 3Kz0.2 simulation yields a smaller maximum q g value at t = 190 min compared to the control simulation.    Interaction with the idealized city also has an impact on the near-surface winds. Figure 19 shows the time-series of the maximum 10-m wind speed in each simulation. In the control simulation, the maximum fluctuates around 10 m s −1 between t = 100 min and t = 170 min. After 170 min, there is a sharp, brief increase to a peak value of 14 m s −1 followed by a slight decrease. This value (14 m s −1 ) represents the smallest peak 10-m wind speed in any of the simulations. In the 0Kz0.5 simulation, the peak maximum is larger compared to the control simulation (24.5 m s −1 vs. 14 m s −1 ) and occurs approximately 10 min earlier. While this peak value is short-lived, values of roughly 15 m s −1 persist until the end of the simulation (not shown). The remaining ∆T = 0 K simulations also exhibit steady maximum winds of roughly 15 m s −1 after interaction with the city, although none have the drastic peak seen in the 0Kz0.5 simulation. Increasing ∆T also results in increased maximum values. For example, in the 3Kz0.2 simulation, the maximum wind speed of 16 m s −1 occurs 36 min prior to the maximum in the control simulation, when the squall line begins interacting with the idealized city. Further increasing ∆T to 5 K results in even stronger 10-m winds. A value of 16 m s −1 occurs during the interaction with the idealized city, and another peak of 21 m s −1 occurs around t = 200 min. In the 7Kz0.2 simulation, the maximum 10-m wind is larger than in the control and 3Kz0.2, but weaker than in 5Kz0.

Discussion
Results indicate that increasing both the magnitude of ∆T and the horizontal gradient in z 0 leads to a deeper and stronger heat island. By many metrics (e.g., maximum updraft velocity, buoyant acceleration, updraft mass flux, graupel mass, 10-m winds), the squall line produced in the control simulation is weaker compared to simulations in which the squall line interacts with a well-defined idealized city. In many instances, increasing the heat island leads to a more substantial modification of the squall line as it passes over the city. However, increasing ∆T or the horizontal z 0 gradient by a defined increment does not always produce an increase in the response (for example, updraft intensity) by a predictable amount, or at all. There are several potential explanations for this finding.
All simulations are identical during the first 80 min of simulation. After this time, small differences begin to emerge. However, the storms do not interact with the city until approximately 70 min later. The most likely cause of these early differences is that some of the inflow air that enters the storm first passes over the idealized city. Analysis of the forward-integrated parcel trajectories reveals that at t = 60 min, when all simulations produce identical storms, air that enters the updraft originates to the left of the idealized cities (not shown). By t = 100 min, when differences between the simulations have emerged, a portion of the parcels originate over and to the right of the idealized cities ( Figure 20). As heat islands of different magnitude develop across the range of simulations, the parcels that pass over the idealized cities will have slightly different temperature, humidity, and wind speed values between the simulations. It is likely that this thermodynamic variability leads to differences in fields such as vertical velocity (see Figure 3) and hydrometeor growth rates (see Figures 16 and 17) well before the storms pass over the idealized city. Even during quasi-steady periods, the updraft experiences oscillations due to the continual generation and decay of new cells at the leading edge of the cold pool (see Figure 3). Fovell and Ogura [42] showed that the frequency of this oscillation can be influenced by microphysical processes. Thus, it is possible that thermodynamic variability within the inflow parcels alters the frequency of the intensification/decay cycle between the simulations. This shift in frequency could impact interpretation of the results when analyzing simulations at a common time and may explain why stronger heat islands do not always indicate stronger storm modification in the analysis presented. parcels that pass over the idealized cities will have slightly different temperature, humidity, and wind speed values between the simulations. It is likely that this thermodynamic variability leads to differences in fields such as vertical velocity (see Figure 3) and hydrometeor growth rates (see Figures  16 and 17) well before the storms pass over the idealized city. Even during quasi-steady periods, the updraft experiences oscillations due to the continual generation and decay of new cells at the leading edge of the cold pool (see Figure 3). Fovell and Ogura [42] showed that the frequency of this oscillation can be influenced by microphysical processes. Thus, it is possible that thermodynamic variability within the inflow parcels alters the frequency of the intensification/decay cycle between the simulations. This shift in frequency could impact interpretation of the results when analyzing simulations at a common time and may explain why stronger heat islands do not always indicate stronger storm modification in the analysis presented.  A second possibility is related to 'chaos seeding' as discussed by Ancell et al. [43]. It is possible that changes to the idealized city generate fast-moving high-frequency waves induced by a numerical roundoff error. These waves result in alterations to the storm structure and lead to non-physical results. While the presence of these waves cannot be completely dismissed, there are a few reasons why they are unlikely to have a significant impact on the results. First, the idealized simulations presented here are of a much shorter duration than those discussed by Ancell et al. [43]. They note that significant changes can be induced "within a day" whereas in the current suite of simulations, the focus is on the structure of the storms after only 2-3 h of integration time. Second, as discussed above, differences between the simulations begin to arise only after the storms ingest air parcels that originate from regions in the domain where the perturbations were placed. Finally, Ancell et al. [43] list several techniques to determine the physical validity of results from perturbation experiments. One such suggestion is to perform a regression analysis of input/response features from an ensemble of simulations. In the context of the current work, it was shown that increasing both the magnitude of the skin temperature perturbation as well as the horizontal z 0 gradient increases the depth and magnitude of the urban heat island. In general, the simulations with the stronger heat islands tended to yield the strongest response in terms of the parameters analyzed (updraft strength, graupel mixing ratio, upward mass flux, 10-m wind speed, etc.), but not always. For example, Figure 13, displays the distribution of buoyant acceleration along parcels entering the storm updraft after interaction with the idealized city. In general, increasing the surface z 0 gradient (while leaving the surface temperature unchanged) or increasing the surface temperature (while leaving the surface z 0 gradient unchanged) leads to an increase in the occurrence of larger buoyant accelerations. Two exceptions to this trend are seen when comparing the 3Kz0.2 simulation to the 3Kz0.5 simulation, and the 3Kz0.2 simulation to the 5Kz0.2 simulation.
An additional exception can be seen in the three simulations with the strongest heat islands (5Kz2.0, 7Kz1.5, and 7Kz2.0). These simulations tend to exhibit a smaller degree of storm modification compared to other simulations with weaker heat islands. In these simulations, deep convection was initiated ahead of the main squall line on the downwind side of the idealized city (e.g., Figures 5,6 and 11). It is possible that the interaction with this new convection weakened the updraft of the main squall line. Figures 10 and 11 indicate that a portion of the air entering the updraft of the main line originated in what should be the cold pool of the new developing convection. Inflow air that travels through these cold pools may have reduced temperature (and potential buoyancy), resulting in slightly weaker updrafts, reduced graupel mass, and weaker 10-m winds than would have occurred if the city-initiated convection did not exist. An additional 7Kz2.0 simulation (not shown) was completed using the same thermodynamic profile, except it was modified so that it contained a capping inversion [44]. In this simulation, new convection was not initiated downwind of the idealized city. Additional work is needed to determine if the results from this current study are valid across a range of convective environments.

Conclusions
The purpose of this work was to simulate city-storm interactions and identify the city characteristics most responsible for structural changes within the storm. To achieve this goal, a suite of idealized two-dimensional simulations was performed in which a simplified representation of a city was inserted into the center of the model domain. In each simulation, the surface temperature and/or surface roughness length within the idealized city were systematically increased relative to the surrounding area by some defined value. While the initial goal was to isolate these two effects, the results show that, in many ways, these two processes work together to decrease low-level stability and promote more vigorous storm updrafts. For example, while increasing the surface temperature over the city does produce a signature similar to an urban heat island, the magnitude and depth of the heat island signature increases with increasing z 0 within the city. In other words, not only does the z 0 gradient have a direct impact on winds, it also increases vertical mixing over the city to create a stronger and deeper urban heat island.
The simulations with the largest heat islands-those with the largest surface temperature and largest roughness gradient-had larger surface-based CAPE and smaller CIN over the idealized city. Compared to the control simulation, the storms in the simulations with heat islands tended to strengthen after they passed over the idealized city. Not only did updrafts become more vigorous, but they became larger, as evident by the increase in cloud top vertical perturbation pressure gradient force as well as the substantial increase in updraft mass flux-an increase that was considerably larger than the observed increase in maximum vertical velocity between simulations.
Interaction with the idealized city also had a strong impact on microphysical processes and near-surface winds. While the storm in the control simulation did experience an increase in graupel mass, the storms in the heat island simulations exhibited a much more substantial increase in graupel during and after interaction with the city. In many cases, modification to the storm structure was not confined to the time period of interaction with the idealized city. Increases in graupel/hail mass, and 10-m wind speeds persisted for more than 30 min after the storms passed over the city. However, the results of the parameter space experiments were neither linear nor monotonic. Increasing the magnitude of the heat island feature did not always lead to a more substantial storm modification. This likely indicates the presence of nonlinear feedback between the environment and in-storm processes. Future work will focus on exploring these feedbacks as well as expanding the simulations to three dimensions.