Modeling the Effect of Desert Urbanization on Local Climate and Natural Dust Generation: A Case Study for Erbil, Iraq

This study uses a suite of meteorological and land-surface models to quantify the changes in local climate and surface dust fluxes associated with desert urbanization. Formulas connecting friction velocity and soil moisture to dust generation are used to quantify surface fluxes for natural wind-blown dust. The models are used to conduct a series of simulations for the desert city of Erbil across a period of rapid urbanization. The results show significant nighttime warming and weak but robust daytime cooling associated with desert urbanization. A slight reduction in near-surface wind speed is also found in the areas undergoing urbanization. These findings are consistent with previous empirical and modeling studies on other desert cities. Numerical models and empirical formulas are used to produce climatological maps of surface dust fluxes as a function of season, and for the preand post-urbanization eras. This framework can potentially be used to bridge the gap in observation on the trends in local dust generation associated with land-use changes and urban expansions.


Introduction
The rapid urban expansion that occurred in desert and semi-arid regions in the past decades has profound impacts on the environment and sustainability [1]. Those dry regions, while already under high environmental stress, have also seen an acceleration of population growth ( [1,2]) that compounds the difficulty of mitigating the environmental impacts. For scientists and stakeholders, it is a major challenge to quantify the changes in critical environmental variables associated with desert urbanization.
Land-use changes associated with desert urbanization affect local weather, water resources, even air quality. Some of those effects are specific to desert cities. For example, while many non-desert cities experience warming in both daytime and nighttime, a weak but robust daytime cooling is known for desert cities from observation ( [3][4][5][6][7]) and numerical modeling ( [8,9]). The arid land surface surrounding a desert city is an ideal ground for the generation of wind-blown dust. The dust fluxes vary with soil moisture and the strength of near-surface turbulence, both of which are affected by urbanization. This aspect has received less attention compared to the effect of temperature. This study aims to use a system of numerical models to quantify the changes in temperature, wind, and natural dust generation associated with desert urbanization.
While long-term observations for surface temperature are available for a small number of desert cities (e.g., [3,4]), similar observations of atmospheric wind and surface dust fluxes are generally lacking for determining the changes in these quantities due to urbanization. Given the limited availability of observations, numerical simulation has been increasingly adopted as an alternative approach for quantifying the effects of urbanization beyond warming and cooling. For example, the changes in wind and diurnal circulation due to urbanization were simulated by [10] for a non-desert city, and [8] for a desert city. These studies did not quantify dust generation. This work adopts the framework of numerical modeling, close to the setting of [8], but further expands it by connecting the meteorological and soil variables to the parameters that determine the surface dust fluxes. The expanded modeling system combines a meteorological model and a land surface model to simulate local wind, temperature, precipitation, and the evolution of soil moisture. These quantities are used as inputs to a set of empirical formulas to quantify natural dust fluxes for pre-and post-urbanization eras. While the approach of combining similar empirical formulas to a meteorological model for the evaluation of dust fluxes has been used before, previous studies (e.g., [11,12]) used it only for short-term weather and air quality forecasts, while this study adopts it for long simulations under land-use changes in the context of desert urbanization.
To test the modeling system, the suite of numerical models is applied to the city of Erbil in northern Iraq. It is a relatively isolated desert city that has undergone a rapid expansion in the turn of the 21st century ( [7,13]). The city was chosen for the first test of the modeling system because of its relatively simple landscape; its urban and suburban areas include only a few types of land cover, and the surrounding region has uncomplicated terrains. Moreover, the city has relatively few high-rise structures but is dominated by low-rise buildings of uniform impervious material ( [13,14]). This allows a relatively simple setting of urban land cover in the numerical models, which allows a clearer analysis and interpretation of the model output. By comparing the numerical simulations for Erbil under preand post-urbanization conditions, we determine the local climate change and the related change in surface dust fluxes induced by urbanization. The relations between local meteorological conditions and dust fluxes are also analyzed.
Although we focus on the simulation for only one city, the modeling framework developed in this study can be applied to other desert cities. The outlook for applications and future generalizations of the modeling system will also be discussed.

Meteorological Model and Boundary Conditions
The numerical simulations in this work use a suite of models: (i) a meso-scale meteorological model that is used to simulate local wind, temperature, and precipitation; (ii) a land surface model that computes the interaction between the land surface and the atmosphere, and the evolution of soil moisture; (iii) a turbulence model that conveys the information of large-scale meteorological fields into the intensity of near surface turbulence; (iv) an empirical model that relates turbulence parameters and soil moisture to dust fluxes. In general, airborne dust can be produced by natural and anthropogenic processes. Due to a lack of pollution inventory for industrial dust generation over desert cities, this study considers only natural wind-blown dust.
Meteorological modeling at urban scale was carried out using the Weather Research and Forecasting (WRF) model [15], Version 3.3.1. All other components of the suite of models are either embedded in WRF or are implemented as extensions to it.
We used nesting of grids in WRF to allow a higher spatial resolution over the target urban area, which is placed in the innermost domain of the model. This domain is nested within larger domains with successively coarser resolutions. To prevent climate drift in seasonal simulations, the lateral boundary conditions for the outermost domain are constrained by large-scale observations from global analysis, specifically the National Centers for Environmental Prediction (NCEP) Final (FNL) data. This is done broadly in the context of dynamical downscaling ( [16][17][18][19][20]).
To process the interaction among the atmosphere, land surface, and sub-surface soil layers, the Noah land surface model [21] was activated within WRF. It calculates the surface fluxes of heat, moisture, and momentum, and simulates the evolution of soil moisture according to a diffusion process and the source and sink due to precipitation and evaporation. The soil moisture at surface will be a key parameter for later calculations of dust fluxes.
Further embedded in the land-surface model is an Urban Canopy Model (UCM), which was also activated in our simulations to modify the surface energy balance and exchanges of fluxes over urban areas. We choose the single-layer UCM developed by [22,23]. The UCM allows users to set urban parameters, such as building heights, roughness length above urban canyons (for momentum equation), urban fraction, and heat capacity of the material in the built environment. Since the details of these parameters over the history of urbanization are not easily obtained, we use a relatively generic setting similar to that in [8] and [9] for desert cities. For example, building height, urban fraction, and heat capacity of roofs were chosen to be 7.5 m, 0.9, and 1 × 10 6 J/m 3 K, respectively. The effects of more complicated spatial arrangements of urban buildings on temperature (e.g., [24,25]) are not considered.
The dynamical core of WRF produces meso-scale meteorological fields at the model grid, which is still too coarse in the boundary layer to resolve turbulence. Since the strength of turbulence, in terms of friction velocity u*, at near-surface is relevant to the calculation of dust fluxes, we also activated a turbulence model embedded in WRF to compute eddy momentum and heat fluxes due to subgrid-scale turbulence. For this purpose, we selected the YSU Non-local-K scheme ( [26,27]).

Calculation of Surface Dust Fluxes
For natural wind-blown dust, surface dust fluxes generally depend on the strength of turbulence, soil moisture, and the type of soil at the surface. These relations were previously determined in laboratory experiments. Our approach is to synthesize existing empirical formulas into a form that can be used to connect WRF model output to dust fluxes. This approach is inspired by the previous study in [12], which used the framework to perform short-term predictions of air pollution. Related surveys can be found in [11,12].
Based on the experimental works of [28][29][30], the vertical mass flux (D f ) for dust particles with a radius less than 10µm is expressed as a function of turbulent friction velocity u*, T , for silt and clay soil (1 − R) 0.13 × 10 −13 (u * ) 3 , if u * ≥ u * T , for sandy soil (1) where u * T is a threshold value, given below, and R is a fractional parameter that depends on land-surface types, as given in Table 1 (compiled after [11,12]) for the USGS 24-category land-cover classification used in WRF.
The threshold value of friction velocity in Equation (1) depends on surface roughness and soil moisture. With a wet surface, a greater strength of turbulence is required to blow dust off the surface. Synthesizing empirical relations developed by [30][31][32], we first compute a baseline value for u * T as where Z 0 is the surface roughness length, which depends on the land-surface type, as given in Table 1. Then, the threshold value, u * T , is computed by where w is the volumetric soil moisture, and w c = 0.0014 C 2 + 0.17 C where C is the fractional content of clay in the soil. In our simulations, we simplify the determination of w c by assuming a 50% clay content, which gives a fixed value of w c = 0.08538. The formulas from the laboratory experiments do not provide the detailed size distribution of the dust. Thus, we restrict the calculations to the bulk (total) dust fluxes. Table 1. U.S. Geological Survey land-use land-cover (LULC) categories and the corresponding surface roughness length (Z 0 ) and dust emission reduction factor (R). The data are compiled after [11] and [12].

The Focus Location and Model Domain
We tested the numerical framework by applying it to a desert city that has undergone a rapid expansion in recent decades. Furthermore, to facilitate a clear interpretation of the numerical results, we chose a city that (i) has a clearly defined boundary and is surrounded by relatively simple pre-urbanization land cover, and (ii) is located over smooth topography so as to minimize the effect of complex terrains. Based on those considerations, we selected Erbil in northern Iraq (centered at approximately 36.2 • N and 44.09 • E) as the target city. The city and its surrounding region are placed in the innermost domain of WRF. It is nested within two more sets of coarser grids. The nested domains are shown in Figure 1. The horizontal resolutions of the innermost, intermediate, and outermost domain are 1 km, 5 km, and 25 km, respectively. The linear dimensions of the three domains are approximately 60 km, 300 km, and 1500 km. (The innermost domain consists of 60 × 60 grid boxes, each with the size of 1 km × 1 km.) Twenty-eight terrain-following vertical levels (the "η-levels" in WRF) are used in all three domains, with the pressure at top of the model set to p = 50 mb (5000 Pa). Figure 2 shows the images from Landsat satellite observation of the city of Erbil and its vicinity in 1987 (left) and 2011 (right). The dark gray pixels correspond to urban land cover. This shows that rapid urban expansion occurred between 1987 and 2011. We chose the land-use maps from these two years to define the pre-and post-urbanization eras in WRF simulations.
Urban expansion over Erbil was generally associated with a conversion from arid types of land to a simple urban-type of land, the latter dominated by concrete and low-rise buildings ( [13,14]). This justifies our choice of using a relatively simple urban land type (see Section 3.2) to represent all "urban" grid cells in WRF. Moreover, the relatively low building heights justify our choice of using the setting of 7.5 m as average building height in the Urban Canopy Model.
Urban Sci. 2020, 4, x FOR PEER REVIEW 5 of 17 "urban" grid cells in WRF. Moreover, the relatively low building heights justify our choice of using the setting of 7.5 m as average building height in the Urban Canopy Model.

Land-Use Land-Cover Maps
The default land-use map in the WRF package is generic (i.e., it is not associated with a specific era but was synthesized from multiple sources across different times) and contains very little urbantype land surfaces associated with the real city of Erbil. Figure 3a shows this default land-use map over the innermost domain of WRF, with urban land-cover marked by a dark blue color. Each grid cell in Figure 3a is 1 km × 1 km. Erbil is located at the center of the domain. Using Landsat observations archived at U.S. Geological Survey ( [33]) and municipal maps of Erbil (extracted from [14]) from the pre-and post-urbanization eras, we created realistic land-use maps over Erbil and its

Land-Use Land-Cover Maps
The default land-use map in the WRF package is generic (i.e., it is not associated with a specific era but was synthesized from multiple sources across different times) and contains very little urbantype land surfaces associated with the real city of Erbil. Figure 3a shows this default land-use map over the innermost domain of WRF, with urban land-cover marked by a dark blue color. Each grid cell in Figure 3a is 1 km × 1 km. Erbil is located at the center of the domain. Using Landsat observations archived at U.S. Geological Survey ( [33]) and municipal maps of Erbil (extracted from [14]) from the pre-and post-urbanization eras, we created realistic land-use maps over Erbil and its

Land-Use Land-Cover Maps
The default land-use map in the WRF package is generic (i.e., it is not associated with a specific era but was synthesized from multiple sources across different times) and contains very little urban-type land surfaces associated with the real city of Erbil. Figure 3a shows this default land-use map over the innermost domain of WRF, with urban land-cover marked by a dark blue color. Each grid cell in Figure 3a is 1 km × 1 km. Erbil is located at the center of the domain. Using Landsat observations archived at U.S. Geological Survey ( [33]) and municipal maps of Erbil (extracted from [14]) from the pre-and post-urbanization eras, we created realistic land-use maps over Erbil and its vicinity for "1987" and "2011", shown in Figure 3b,c, which were used to override the default map in WRF. This was applied only to the city and its vicinity. Elsewhere, the default land-use map in WRF was retained. This design is consistent with our idea of running "twin experiments", explained below, with the only difference between the two runs coming from the surface boundary condition related to urbanization. Figure 3b,c show that the urban type of land (dark blue) expanded significantly from 1987 (~94 km 2 ) to 2011 (~284 km 2 ). This detail of land-use change is well resolved by the 1 km × 1 km grid (imposed in the background of Figure 3a-c) in the innermost domain of WRF. A few bright yellow grid boxes that emerge in western Erbil in 2011 are associated with a new park that was constructed between 1987 and 2011. There, the land-surface type is mixed grassland and shrubland instead of urban. The non-urban area to the east of Erbil is dominated by dry cropland (blue in Figure 3) with a greater potential of dust generation in contrast to the cropland/woodland (green) and shrubland (tan) that dominate the non-urban areas to the west of Erbil.
Urban Sci. 2020, 4, x FOR PEER REVIEW 6 of 17 vicinity for "1987" and "2011", shown in Figure 3b,c, which were used to override the default map in WRF. This was applied only to the city and its vicinity. Elsewhere, the default land-use map in WRF was retained. This design is consistent with our idea of running "twin experiments", explained below, with the only difference between the two runs coming from the surface boundary condition related to urbanization. Figure 3b,c show that the urban type of land (dark blue) expanded significantly from 1987 (~94 km 2 ) to 2011 (~284 km 2 ). This detail of land-use change is well resolved by the 1 km × 1 km grid (imposed in the background of Figure 3a-c) in the innermost domain of WRF. A few bright yellow grid boxes that emerge in western Erbil in 2011 are associated with a new park that was constructed between 1987 and 2011. There, the land-surface type is mixed grassland and shrubland instead of urban. The non-urban area to the east of Erbil is dominated by dry cropland (blue in Figure 3) with a greater potential of dust generation in contrast to the cropland/woodland (green) and shrubland (tan) that dominate the non-urban areas to the west of Erbil.

Lateral Boundary Conditions and Dynamical Downscaling
With the model domain and land-surface boundary conditions given in the preceding sections, we conducted a series of numerical simulations in the manner of "twin experiments". To extract the effect of urbanization, a pair of runs were performed with identical external forcing (seasonally and diurnally-varying solar radiation) and boundary conditions, except the land-use map in the surface boundary condition. For this purpose, the lateral boundary conditions were constructed from the large-scale observations for a generic year, chosen as 2000. Hereafter, the two runs named "1987" and "2011" in a twin experiment refer to the simulations with the same generic lateral boundary conditions but with the land-use map in WRF replaced by those from 1987 and 2011, respectively. With this setup, one can attribute the difference between the two runs to the effect of urbanization through the changes in the surface boundary conditions.
Since local climate and dust generation vary with season, we conducted two runs for each of the 1987 and 2011 cases. The "summer" simulation lasted for 30 days in the month of July (from July 1 to 30), and the "winter" simulation for 30 days in the month of January (from January 1 to 30), as summarized in Table 2. In each of the month-long simulations, the lateral boundary conditions for the outermost domain were constrained by four daily observations from global analysis, using the FNL data set. The output of major meteorological variables (e.g., temperature and precipitation), friction velocity produced by the turbulence model, and soil moisture produced by the land-surface model were all saved at 3-h intervals. Each run has a total of 240 3-hourly outputs which are used for further analyses.

Cross Validation with Reanalysis
Climate studies over desert cities have long been hindered by a lack of reliable observational networks, particularly at the sub-meso and urban scales. This motivated us to develop a numerical modeling framework as a complementary tool for research. Nevertheless, we will briefly compare a basic simulation for the Erbil region with reanalysis, which is treated as the best available (if heavily interpolated) representation of the raw observation. We focus on thermodynamic variables, as the velocity field generally has strong spatial and temporal variations at the sub-mesoscale which would not be properly represented in reanalysis. Using Climate Forecast System Reanalysis (CFSR) data [34] as the basis, we selected two grid locations closest to Erbil: L1 at (36.06 • N and 43.75 • E), and L2 at (36.06 • N and 44.06 • E). L2 is approximately 10 km south of Erbil; L1 is at the same latitude as L2 but 25 km west of it. Since the resolution of the WRF model is higher that of the reanalysis, for the comparison to the latter, we took the average of the model output over multiple WRF grid points surrounding L1 and L2, respectively. Figure 4 compares the daily values from the CFSR reanalysis (red) and the WRF simulation (black) for L1 (left column) and L2 (right column) over the month of January 2000. In each column, top to bottom, are 2 m air temperature, relative humidity, daily rainfall, and 10 m wind speed. For temperature and humidity, the WRF simulation is consistent with CFSR in the temporal variation. For precipitation, strong events are present in both WRF and CFSR, while discrepancies are seen for some weak events. Note that precipitation in the reanalysis is not directly constrained by observation but is generated by the global model which depends on physical parameterization. A more notable systematic difference between WRF and CFSR is in the 10 m wind speed, with CFSR having a lower magnitude. Nevertheless, consistency is still found in the temporal variation (timing of occurrences of maximum wind speed) between the two. Like precipitation, velocity in the surface boundary layer is not directly constrained by observation in reanalysis. Given that the CFSR global model has a lower resolution than in the WRF model, it is not surprising that events with maximum wind gusts are smoothed out and more muted in the former.
The overall consistency between the time series from WRF and CFSR reanalysis affirms the soundness of our setup of time-varying boundary conditions for dynamical downscaling. To more accurately assess the differences in the details for precipitation and surface wind speed, we would require local observational data at the sub-meso and urban scale, which we hope will become available in the future.

Effect of Urbanization on Meteorological Fields
Previous numerical studies on the effect of desert urbanization have focused on surface air temperature. The well-known "urban heat island" effect was shown to generally hold for desert cities in nighttime, although this is accompanied by weak daytime cooling ( [8,9]).
To compare with previous studies, Figure 5 shows the 3-hourly time series of 2 m air temperature, averaged over the grid boxes where urbanization occurred between 1987 and 2011, for the 1987 and 2011 simulations. Figure 5a,b are for January and July, respectively, and the runs with the 1987 and 2011 land-cover are colored in black and red, respectively. In both figures, even on a day-to-day basis, one observes a systematic increase in nighttime temperature from 1987 to 2011 while the temperature difference is smaller in daytime. (Nevertheless, daytime cooling can be identified more clearly after monthly average, as shown below). This is consistent with the findings from previous studies for different desert cities (e.g., [3,4,8,9]). Figure 6a,b show the maps of the difference, defined as "2011 minus 1987", in nighttime 2 m air temperature averaged over the entire month of January and July, respectively. Warming over the newly established urban area is readily identified in both plots, with a magnitude of warming at 1.5 • C and 2 • C for winter and summer, respectively. (Here, we use the average of midnight and 3 a.m. local time to define "nighttime", as the maximum of nighttime warming usually occurs in early morning). Figure 6c,d are similar to Figure 6a,b but for the temperature difference in daytime (defined as the average of noon and 3 p.m.). A weak but robust signal of daytime cooling (0.25 • C and 0.45 • C for winter and summer, respectively) due to urbanization is found over the newly established urban areas. The pattern of "strong warming in nighttime and weak cooling in daytime" occurs in both winter and summer, but with the signal stronger in summer. Overall, these findings are consistent with previous numerical studies on desert urbanization. Furthermore, note that anthropogenic effects, such as increased irrigation within the city, or anthropogenic production of heat, are not explicitly included in our study. As explained by [8,9], daytime cooling in the type of simulations we performed arises mainly from an increased effective surface area over the city (due to presence of buildings; equivalent to the "shadow effect").
A relatively small number of studies have examined the effect of urbanization on the near-surface winds. A previous study [8] found a reduction in the strength of diurnal circulation over Las Vegas (a desert city) due to urbanization, while another study [35] found a similar reduction in low-level wind speed over a tropical non-desert city of Muar, Malaysia. The changes in wind field could, in turn, affect the temperature by altering the "ventilation effect" [8]. To compare with existing works, Figure 7 shows the 3-hourly time series of the 10 m wind speed averaged over the area that underwent urbanization, in the same format as Figure 5. Figure 7a,b are for winter and summer, respectively, and black and red curves are for the 1987 and 2011 cases, respectively. From the figure, one observes a slight downward shift of wind speed from 1987 to 2011. This reduction in wind speed can be more clearly identified in the maps of the monthly mean, as shown in Figure 8 (in the same format as Figure 6). It is found that urbanization leads to a reduction in surface wind speed in both winter and summer, and both daytime and nighttime. This is similar to the findings of the numerical study by [8] for Las Vegas.

Surface Dust Fluxes
Using the empirical formulas in Section 2.2, the 3-hourly outputs of friction velocity u* and soil moisture w from the WRF simulations were used to compute the surface dust fluxes over the innermost model domain. Figure 9a,b show the maps of monthly cumulative dust fluxes from the 1987 case for winter and summer, respectively. The outline of the border of urban areas of Erbil in 1987 is superimposed to the maps. Figure 9c,d are similar to Figure 9a,b but for the 2011 case.
As expected, the conversion from arid-type natural land to urban-type land leads to suppression of dust generation; in Figure 9, the areas with minimum dust fluxes are those covered by urban land. Nevertheless, after urban expansion, more people are now living closer to the "hot spots" just outside

Surface Dust Fluxes
Using the empirical formulas in Section 2.2, the 3-hourly outputs of friction velocity u* and soil moisture w from the WRF simulations were used to compute the surface dust fluxes over the innermost model domain. Figure 9a,b show the maps of monthly cumulative dust fluxes from the 1987 case for winter and summer, respectively. The outline of the border of urban areas of Erbil in 1987 is superimposed to the maps. Figure 9c,d are similar to Figure 9a,b but for the 2011 case.
As expected, the conversion from arid-type natural land to urban-type land leads to suppression of dust generation; in Figure 9, the areas with minimum dust fluxes are those covered by urban land. Nevertheless, after urban expansion, more people are now living closer to the "hot spots" just outside

Surface Dust Fluxes
Using the empirical formulas in Section 2.2, the 3-hourly outputs of friction velocity u* and soil moisture w from the WRF simulations were used to compute the surface dust fluxes over the innermost model domain. Figure 9a,b show the maps of monthly cumulative dust fluxes from the 1987 case for winter and summer, respectively. The outline of the border of urban areas of Erbil in 1987 is superimposed to the maps. Figure 9c,d are similar to Figure 9a,b but for the 2011 case.
As expected, the conversion from arid-type natural land to urban-type land leads to suppression of dust generation; in Figure 9, the areas with minimum dust fluxes are those covered by urban land.
Nevertheless, after urban expansion, more people are now living closer to the "hot spots" just outside the city limit where strong dust generation occurs. These hot spots for dust generation are located mainly to the west of the city, where the dominant land surface type is dry cropland (see Figure 3b,c). To the east of the city, where land cover is dominated by cropland/woodland and shrubland, dust fluxes are relatively weak. This is due to those two land surface types having less clay content (thus, greater values of the reduction parameter R in Table 1; both have R = 0.7, compared to R = 0.4 for dry cropland).
Another reason for the stronger dust production to the east of Erbil is that a condition of stronger turbulent generally persists in that region. Figure 10 shows the monthly mean of friction velocity u * from WRF simulations for the 1987 case, for January and July. Both seasons exhibit a clear contrast in friction velocity which is stronger in the east and weaker in the west. Over the hot spots for dust generation, the dust fluxes have greater values in summer. This is due to both a higher friction velocity and drier land surface in summer. Figure 11 shows the difference in the cumulative dust fluxes between 1987 and 2011 for winter and summer. For this particular set of maps, since dust generation is suppressed over the newly established urban areas in 2011, we used the reverse definition of "1987 minus 2011" to define the difference. (Thus, a positive value means more dust generation in 1987). As expected, the areas with positive differences in dust generation approximately coincide with those that underwent urbanization between 1987 and 2011. Nevertheless, even within those areas, there are nontrivial spatial variations in the magnitude of the differences. This illustrates the usefulness of the meteorological model simulations, as the spatial variations of meteorological fields help create these detailed patterns.
Urban Sci. 2020, 4, x FOR PEER REVIEW 12 of 17 the city limit where strong dust generation occurs. These hot spots for dust generation are located mainly to the west of the city, where the dominant land surface type is dry cropland (see Figure 3b,c).
To the east of the city, where land cover is dominated by cropland/woodland and shrubland, dust fluxes are relatively weak. This is due to those two land surface types having less clay content (thus, greater values of the reduction parameter R in Table 1; both have R = 0.7, compared to R = 0.4 for dry cropland). Another reason for the stronger dust production to the east of Erbil is that a condition of stronger turbulent generally persists in that region. Figure 10 shows the monthly mean of friction velocity u * from WRF simulations for the 1987 case, for January and July. Both seasons exhibit a clear contrast in friction velocity which is stronger in the east and weaker in the west. Over the hot spots for dust generation, the dust fluxes have greater values in summer. This is due to both a higher friction velocity and drier land surface in summer. Figure 11 shows the difference in the cumulative dust fluxes between 1987 and 2011 for winter and summer. For this particular set of maps, since dust generation is suppressed over the newly established urban areas in 2011, we used the reverse definition of "1987 minus 2011" to define the difference. (Thus, a positive value means more dust generation in 1987). As expected, the areas with positive differences in dust generation approximately coincide with those that underwent urbanization between 1987 and 2011. Nevertheless, even within those areas, there are nontrivial spatial variations in the magnitude of the differences. This illustrates the usefulness of the meteorological model simulations, as the spatial variations of meteorological fields help create these detailed patterns.

Influence of Meteorological Conditions on Dust Generation
The empirical formulas for dust generation given in Section 2.2 quantify the impact of meteorological conditions on the dust fluxes. Qualitatively, over the same type of land surface, dust fluxes increase under stronger turbulence and lower soil moisture. To illustrate how the two factors influence the dust fluxes under arid climate over Erbil, we constructed time series of relevant variables at a selected WRF grid box located at approximately (44.09° E, 36.17° N), which is at the edge of Erbil where the land turned from non-urban in 1987 to urban in. Figure 12a shows the 3hourly time series, from the January 1987 case, of friction velocity (u*, black) along with the threshold value ( * , red) for dust production. Figure 12b shows the time series of soil moisture. The few upticks of soil moisture correspond to events of precipitation. Figure 12c shows the local dust fluxes (Df).
From Figure 12, the local u* frequently exceeds the threshold value that triggers dust production. While the threshold of u* increases with soil moisture, this effect is minor, as the red curve in Figure  12a only goes up slightly with the upticks in soil moisture. This is due to the fact that the surface over the desert is generally very dry. Even when it rains, soil moisture only goes up slightly and quickly drops back to the baseline dry value. As such, dust generation is largely controlled by the friction

Influence of Meteorological Conditions on Dust Generation
The empirical formulas for dust generation given in Section 2.2 quantify the impact of meteorological conditions on the dust fluxes. Qualitatively, over the same type of land surface, dust fluxes increase under stronger turbulence and lower soil moisture. To illustrate how the two factors influence the dust fluxes under arid climate over Erbil, we constructed time series of relevant variables at a selected WRF grid box located at approximately (44.09° E, 36.17° N), which is at the edge of Erbil where the land turned from non-urban in 1987 to urban in. Figure 12a shows the 3hourly time series, from the January 1987 case, of friction velocity (u*, black) along with the threshold value ( * , red) for dust production. Figure 12b shows the time series of soil moisture. The few upticks of soil moisture correspond to events of precipitation. Figure 12c shows the local dust fluxes (Df).
From Figure 12, the local u* frequently exceeds the threshold value that triggers dust production. While the threshold of u* increases with soil moisture, this effect is minor, as the red curve in Figure  12a only goes up slightly with the upticks in soil moisture. This is due to the fact that the surface over the desert is generally very dry. Even when it rains, soil moisture only goes up slightly and quickly drops back to the baseline dry value. As such, dust generation is largely controlled by the friction Figure 11. The difference in the monthly cumulative dust flux, defined as "1987 case minus 2011 case". Positive means more dust is produced in the 1987 case. (a) Winter (January); (b) summer (July).

Influence of Meteorological Conditions on Dust Generation
The empirical formulas for dust generation given in Section 2.2 quantify the impact of meteorological conditions on the dust fluxes. Qualitatively, over the same type of land surface, dust fluxes increase under stronger turbulence and lower soil moisture. To illustrate how the two factors influence the dust fluxes under arid climate over Erbil, we constructed time series of relevant variables at a selected WRF grid box located at approximately (44.09 • E, 36.17 • N), which is at the edge of Erbil where the land turned from non-urban in 1987 to urban in. Figure 12a shows the 3-hourly time series, from the January 1987 case, of friction velocity (u*, black) along with the threshold value (u * T , red) for dust production. Figure 12b shows the time series of soil moisture. The few upticks of soil moisture correspond to events of precipitation. Figure 12c shows the local dust fluxes (D f ).
From Figure 12, the local u* frequently exceeds the threshold value that triggers dust production. While the threshold of u* increases with soil moisture, this effect is minor, as the red curve in Figure 12a only goes up slightly with the upticks in soil moisture. This is due to the fact that the surface over the desert is generally very dry. Even when it rains, soil moisture only goes up slightly and quickly drops back to the baseline dry value. As such, dust generation is largely controlled by the friction velocity.
The story is similar in summer (not shown). In fact, at the same location, it hardly rains throughout the month of July, and the time series of soil moisture and threshold value of u* would remain flat all month. The friction velocity is generally greater in summer, corresponding to a greater value of the dust flux.
Urban Sci. 2020, 4, x FOR PEER REVIEW 14 of 17 velocity. The story is similar in summer (not shown). In fact, at the same location, it hardly rains throughout the month of July, and the time series of soil moisture and threshold value of u* would remain flat all month. The friction velocity is generally greater in summer, corresponding to a greater value of the dust flux.

Concluding Remarks
This study uses a suite of meteorological and land-surface models to quantify the effects of urbanization on local climate and surface dust fluxes for an isolated desert city. From a set of simulations for Erbil, Iraq, a general pattern of significant nighttime warming and weak but robust daytime cooling associated with urbanization was found in 2 m air temperature. Urbanization also leads to a slight reduction of near-surface wind speed at both daytime and nighttime. These effects exist in both winter and summer, with the temperature effect stronger in summer. These conclusions are consistent with, and further strengthen, the results from previous numerical modeling studies on desert urbanization.
The new aspect of modeling accomplished in this work is a further connection of meteorological and soil variables to a set of empirical formulas for quantifying surface dust fluxes. Using friction velocity and soil moisture from the meteorological model as the main input parameters, the empirical formulas produce climatological maps of surface dust fluxes as a function of season, for the pre-and post-urbanization eras. Since field observations of dust fluxes are not available for Erbil, the results from the numerical simulations serve as the virtual observations from which we gain the first estimate of the effect of urbanization on dust, and further insights into the interaction between

Concluding Remarks
This study uses a suite of meteorological and land-surface models to quantify the effects of urbanization on local climate and surface dust fluxes for an isolated desert city. From a set of simulations for Erbil, Iraq, a general pattern of significant nighttime warming and weak but robust daytime cooling associated with urbanization was found in 2 m air temperature. Urbanization also leads to a slight reduction of near-surface wind speed at both daytime and nighttime. These effects exist in both winter and summer, with the temperature effect stronger in summer. These conclusions are consistent with, and further strengthen, the results from previous numerical modeling studies on desert urbanization.
The new aspect of modeling accomplished in this work is a further connection of meteorological and soil variables to a set of empirical formulas for quantifying surface dust fluxes. Using friction velocity and soil moisture from the meteorological model as the main input parameters, the empirical formulas produce climatological maps of surface dust fluxes as a function of season, for the pre-and post-urbanization eras. Since field observations of dust fluxes are not available for Erbil, the results from the numerical simulations serve as the virtual observations from which we gain the first estimate of the effect of urbanization on dust, and further insights into the interaction between meteorological conditions and dust generation over the desert city. In particular, for the case of Erbil, it was found that the dust fluxes are controlled mainly by the strength of near-surface turbulence, while soil moisture plays a relatively minor role.
The modeling system tested in this work is general enough that it can be readily adopted for quantifying the effect of urbanization for other cities. The required inputs from local stakeholders for a specific city are the land-use maps in the surface boundary condition and the urban parameters (building height, fraction of urban-type surface, etc.) in the urban canopy model. The modeling system can be used to evaluate not only the difference between past and present climate states, but also the difference between present and future given a projection of future land-use changes based on city planning.
For future improvements of the system, note that the existing empirical formulas for computing dust fluxes only return the bulk mass flux for all particles with radius <10 µm. If refined formulas with additional information on dust size distribution become available from updated laboratory experiments, such formulas can be readily inserted into our framework as an upgrade. Knowing the spectrum of dust particles will be useful for computing the transport of dust by the local atmospheric circulation, another component that can be added to the system in future developments.