Modelling atmospheric dry deposition in urban areas using an urban canopy approach

This section doesn’t provide the model-observation comparison. In fact, it summarized the range and uncertainty of measurements in previous studies. I suggest shortening this section to one paragraph and insert it into Section 5.2. A: We agree that the title of this section is not really consistent with its content. We think however that this content is valuable. Therefore it has been moved to an appendix titled “Overview of 10


Introduction
Although urban areas currently occupy only a few percent of the Earth's surface (2.8 % in 2011(2.8 % in , Martine, 2011, more than half of the world's population lives in urban areas. This figure reaches at least 80 % in Europe, North America and Japan (Elvidge et al., 2004, as cited by Oleson et al., 2008Oleson et al., , p. 1039) and urban areas are expected to increase in the future (Shepherd, 2005). Consequently, the health and environmental impacts of pollutants within these urban areas are of great concern in air quality studies. The deposition fluxes of air pollutants have rarely been modelled within urban areas.
Historically, atmospheric deposition studies have focused mostly on remote areas to assess the potential impacts on ecosystems of acid deposition and nitrogen loading, or the potential impacts on human health of pollutants such as mercury or persistent organic pollutants, which bioaccumulate in the food chain (e.g. 51 % of cereals consumed in France in 2004 contained pesticides, de Jaeger et al., 2012). Therefore, current atmospheric deposition models may not be suitable to simulate deposition fluxes in urban areas, which include complex surface geometries and diverse land use types. Atmospheric deposition in urban areas is a topic of current interest for several reasons. For example, there is a growing interest for urban horticulture (Säumel et al., 2012) and green roofs (Yang et al., 2008), and vegetation may be adversely affected by atmospheric pollutant deposition. Air pollutant deposition on buildings and other surfaces may lead to soiling and degradation of their surfaces, thereby leading to cleaning or replacement costs as well as loss of architectural/cultural value. Furthermore, atmospheric deposition contributes to the contamination of storm water and the mobilisation of pollutants by water runoff depends on the surface type and configuration. Both wet and dry processes contribute to atmospheric deposition. Models of wet deposition do not depend on the surface type, and can, therefore, apply to all types of areas, including urban areas. On the other hand, dry deposition depends strongly on the surface type and there is a need to develop dry deposition models that take into account the characteristics of urban areas (Maro, 2014;Jonsson et al., 2008). Currently, dry deposition models are too simple for application to the urban environment. Their classical approaches (Wesely and Hicks, 2000;Petroff et al., 2008a), which are inherited from semi-empirical models, were developed for deposition over vegetated surfaces, bare soil or water, and therefore they fail to represent the complexity of the dry deposition processes over an urban canopy. We present here the development and initial application of a dry deposition model for the urban environment.

Brief historical review of the dry deposition velocity
The mass transfer of pollutants between the air and exposed surfaces is controlled by a wide range of chemical, physical and biological processes, which may interact among each others. However, for the sake of simplicity, the concept of deposition velocity was introduced. Gregory (1945) first introduced this concept as the ratio of the flux F of an air pollutant towards a surface measured at a reference height z ref and its concentration c measured at the same height, leading to the following formulation This formulation allows one, through the knowledge of v d , to estimate the dry deposition flux F from the airborne concentration c in a three-dimensional air quality model: where x, y are the horizontal coordinates and t the time. For gases, the dry deposition velocity is generally computed from a formulation analogous to Ohm's law in electrical circuits (Wesely and Hicks, 2000), e.g.
where R a , R b and R c are resistances to mass transfer. Each resistance represents the process that predominantly governs mass transfer from the air towards the surface. For the turbulent regime of the surface layer, the aerodynamic resistance, R a , represents the resistance to turbulent mass transfer. It has the same value for all substances and depends solely on the atmospheric flow. R b represents the quasilaminar resistance to mass transfer via molecular diffusion through the thin laminar layer of air (a few millimetres) just above the surface. R c is called the surface resistance; it takes into account the interaction processes (adsorption, absorption, chemical reaction, etc.) between the surface and the substances being deposited. R b and R c depend on the substance characteristics. For particles, the latter two resistances, R b and R c , are generally replaced by a lumped surface resistance, R s , (e.g. Slinn, 1982) and the contribution of gravitational settling that becomes relevant for the coarser particles is integrated within the v d formulation (Venkatram and Pleim, 1999).
In this work, we focus on the aerodynamic resistance, because it depends mainly on the atmospheric flow characteristics, and therefore is strongly influenced by the urban canopy.

Existing urban canopy models
Numerous urban modelling schemes have been developed in the past decade (e.g. Brown, 2000) to approximate the effect of the local-scale urban elements on drag, heat flux and the radiative budget. Large-scale numerical models do not have the spatial resolution needed to represent fluid dynamics at the scales relevant to the built urban environment. Several reviews of urban models are available (e.g. Brown, 2000;Masson, 2006;Grimmond et al., 2010Grimmond et al., , 2011. For example, Masson (2006) considers three general categories of urban parametrisations: -Empirical models: these models are based on observations and represent the behaviour of the urban canopy using statistical relations.
-Vegetation models: these models have been modified to fit to urban characteristics.
-Urban canopy models: these relatively recent models include a representation of the urban canopy in the dynamic flow equations.
These latter urban canopy models are based on simple geometries, but are nevertheless appropriate to represent the main aerodynamic and thermal characteristics of urban areas. However, they have so far been intended to parametrise the momentum and heat transfer processes, not dry deposition of atmospheric pollutants.
The combination of previously existing concepts allows us to propose here a novel approach to model dry deposition of atmospheric pollutants in an urban canopy. It is based on the urban canyon concept of Oke (1988). The modelling concept is based on a single infinitely long road, bordered by two facing buildings, which are treated separately. It accounts for local effects of buildings through the use of a local mixing length and key parameters characteristic of the urban canopy. Three different flow regimes are distinguished in the urban canyon according to its height-to-width ratio: isolated roughness flow, wake interference flow and skimming flow (Oke, 1987). The turbulence scheme used in the classical roughness-length approach using the wind velocity to parametrise turbulent motions is modified to make it suitable for the urban canopy. This approach provides spatially distributed dry deposition fields within the urban canopy, which cannot be obtained from the roughness-length model.
We summarise first the formulation of the roughnesslength model. Next, we describe the subgrid model developed here and present the dry deposition flux for the different flow regimes and surface types. Finally, simulations are conducted to compare the dry deposition fluxes obtained with this model and the roughness-length model, as well as to investigate the sensitivity of the model results to several key parameters.
where x i are the Cartesian coordinates, u i and u i are, respectively, the mean and the fluctuating components of the wind velocity in the direction x i , c is the fluctuating component of the concentration c, and D is the molecular (for gases) or Brownian (for particles) diffusivity. S represents other sources or sinks of the pollutant. A closure problem arises because of the non-linear term u i c . The analogous terms in the Reynolds-averaged Navier-Stokes (RANS) equations are known as the Reynolds stress: In order to close the system of RANS equations, Boussinesq introduces the turbulent momentum diffusivity to provide a widely used relationship between the Reynolds stress and the mean terms of the flow fields (e.g. Schmitt, 2007). In the surface layer (at least in the upper part), this hypothesis allows one to formulate the turbulent momentum flux as follows: where ν t is the turbulent momentum diffusivity (eddy viscosity), u and u are respectively the mean and the fluctuating components of the wind velocity parallel to the considered surface, w is the fluctuating component of the normal wind velocity and z is the coordinate along the normal to the surface. By analogy, the first-order closure scheme for mass transfer, also called K-theory, leads to the following formulation of the vertical turbulent mass flux: where K c t is the turbulent mass diffusivity. The only available framework, which allows one to express the deposition velocity as a function of resistances, assumes that the vertical mass flux is constant. This assumption is valid far away from a rough surface. The vertical mass flux is the sum of the turbulent mass flux F c t , which dominates in the atmospheric turbulent layer, and the molecular diffusion mass flux, F c D , which dominates only in the quasi-laminar sublayer near the surface. Thus, when calculating the aerodynamic resistance, F c ≈ F c t . Subsequently, where z b is the height at which turbulent motions stop governing mass transfer compared to Brownian motion. Subsequently, the aerodynamic resistance may be expressed as follows: Although ν t is reasonably well known, this is not the case for K c t . A standard approach consists in relating K c t to ν t through the following ratio: where Sc t is the turbulent Schmidt number.
In the surface layer, well above the canopy, the standard assumption that the eddy diffusivities for concentration and temperature are equal to the turbulent viscosity (i.e. K ϕ t = ν t , where ϕ may be either the concentration or the temperature) is generally accepted (Businger, 1986). Within the roughness sublayer (RSL) (generally defined as the sublayer where the standard flux-gradient relationships fail), these eddy diffusivities are modified for temperature (turbulent Prandtl number, Pr t , different from 1) and concentration (Sc t = 1). Petroff (2005) 1 partly explains the difference between the turbulent transport of momentum and that of scalars (temperature and concentration) by the influence of the canopy on the flow fields, such as the production of Rayleigh instability for the temperature (Raupach et al., 1991).
Within urban areas, Sini et al. (1996) chose Pr t to be equal to 0.7. Concerning Sc t , very few studies have been conducted, especially in urban areas. Tominaga and Stathopoulos (2007) showed that Sc t should be close to 0.3 around a single building. However, they argue that a large number of buildings should produce additional turbulence, which would lead to a greater value of an effective Sc t . Because of the lack of studies, Sc t is generally chosen equal to unity. This assumption impacts all the deposition models considered here in the same way.
The Prandtl mixing-length theory is a widely used model to parametrise the turbulent eddy viscosity in the atmospheric surface layer. It allows one to express the turbulent viscosity as follows: where l m a characteristic mixing length for turbulent motion. It leads to the following aerodynamic resistance formulation, used in most operational air quality models (e.g. Zhang et al., 2001): where z 0 is the roughness length. However, the Prandtl mixing-length theory leads to formulations that are only valid in a region far enough from the surface so that viscous effects can be neglected (e.g. Schlichting and Gersten, 2000). For very rough surfaces (forest, urban areas, etc.), the influence of the surface can be significant at distances that are not negligible (up to several canopy heights, e.g. Thom et al., 1975;Raupach, 1979). This layer is usually known as the RSL. Subsequently, the introduction of a zero-plane displacement height, d, is a commonly used approximation. The resulting formulation is considered satisfactory to represent the dry deposition flux as a sink for atmospheric concentrations. However, this model does not provide any detailed information on the dry deposition processes occurring inside the urban canopy.

Urban canopy model
The model described here is developed for use in threedimensional gridded air quality models and is designed to simulate the transfer of pollutants from the atmosphere to urban surfaces. It is a bulk approach, developed using a subgrid parametrisation. Thus, only the lowest grid layer will be investigated. In air quality models, the lowest model layer is generally between 25 and 50 m high (e.g. van Loon et al., 2007), although heights as low as 14 to 25 m have been reported in recent applications (Solazzo et al., 2013). It is assumed here that the height of the lowest model layer is at least twice that of the urban canopy. The currently available roughness-length models use some urban canopy parameters (roughness length, displacement height, etc.) to estimate dry deposition in urban areas but it is not designed to reproduce the flow fields within the urban canopy. Here, we use the canyon concept developed by Nunez and Oke (1977). The urban canyon consists of a single road, bordered by two facing buildings, which can be treated separately. The individual shapes of individual buildings are not taken into account and only spatially averaged characteristics of the urban area (mean building height h, canyon width W , etc.) are used. Any road orientation is possible and exists with the same probability.
The flow fields depend on the canyon geometry. The range of canyon geometries is split into three different flow regimes depending on the height-to-width ratio of the canyon: -In a very narrow canyon, a vortex can develop within the canopy, leading to a recirculation region (noted as r in the variable subscript), similar to a cavity flow, which is called skimming flow.
-If the canyon is large enough, a second region, the ventilation region (noted as v in the variable subscript), appears downwind of the recirculation region. The flow pattern is called isolated roughness flow.
-Between these two cases, the downwind buildings leads to a ventilation region that does not extend down to the ground. This flow pattern is called wake interference flow.
The boundaries between these two regions still need to be defined. In most models using this approach, the shape of the recirculation region is a trapezoid (e.g. see Fig. 2). According to the review by Harman et al. (2004), measurements show that the maximum length of the recirculation region (the base of the trapezoid, W r ) is proportional to the height of the building, h. Harman et al. (2004) show that the ratio W r h depends on the turbulence level in the boundary layer and the shape of the buildings and roofs. For a cubical array of buildings (a hypothesis assumed by Macdonald et al., 1998, for the calculation of the displacement height d), Castro and Robins (1977) proposed W r h ≈ 2. On the other hand, Oke (1988) Okamoto et al. (1993) described a two-dimensional geometry, which resembles a realistic urban area, and recommended W r h ≈ 3.5. Here, we selected W r = 3 h. The sensitivity of the model to this value is tested in Sect. 5.4.2.
The three flow regimes are then split according to the length of the flow regions (in particular the recirculation region): -For narrow canyons ( Fig. 1), W < W r 2 , i.e. h W > 2 3 , which corresponds to the skimming flow regime.

Parametrisation of turbulence within the urban canopy
As already stated, the standard flux-gradient relationships fail to reproduce the mean flow and concentration profiles within and above an urban canopy. Applying K-theory to the transport of pollutants may be even more problematic than its application to momentum, because the length scales involved in the transport of pollutants are even smaller than those involved in the transport of momentum.
Numerous schemes have been developed for momentum, such as non-local closure schemes (e.g. probability density function theory Pope, 2000, or the transilient theory from Stull, 1984). Concerning pollutant concentrations, Raupach (1989) developed an alternative to K-theory with its localised near-field theory (LNF) within vegetative canopies. This latter theory splits the pollutant transport into two components: advection from near-field sources and diffusion from far-field contributions.
Such approaches are generally considered too demanding in terms of computational requirements and/or input data (e.g. source or sink distribution) for routine application in air quality modelling. Therefore, all these constraints (computational costs, lack of available data, etc.) point out the need for a simple model (such as flux-gradient relationships) to predict dry deposition fluxes above and within the canopy.
This work aims to develop a revised flux-gradient relationship, based on an improved length scale of turbulence compared to that used in the roughness-length model (Sect. 3.2.1), coupled to a realistic representation of the wind speed profile within the canopy (Sect. 3.2.2).

Urban mixing length
First, we improve the characteristics of the mixing length compared to that used in the roughness-length model. The impact of buildings can be taken into account by introducing a new mixing length. The roughness elements, such as buildings, generate turbulent wakes, and the size of resulting eddies is known to be related to the dimensions of these roughness elements.
Following Coceal and Belcher (2004), the general form of the mixing length will be deduced from the following two extreme cases: -If the urban canopy has low building density, turbulence should not be affected significantly by the urban canopy.
In this case, turbulent eddies are blocked mostly by the ground and the mixing length, l m , follows the law of the wall profile: l m = κz, where z is the distance to the surface and κ is the von Kármán constant (taken here to be 0.41).
-If the urban canopy is very dense, the large eddies above the urban canopy break at the top of the canopy. Raupach et al. (1996) show that the dominant eddies within a vegetation canopy are mostly produced from mixinglayer instability of the shear layer, which is created at the top of the canopy. The mixing length in a very dense canopy, l c , is then assumed to be constant with height in order to reflect this behaviour, controlled by the thickness of the shear layer. It is then expected to depend on the mean height of buildings. Coceal and Belcher (2004) proposed to interpolate between these two behaviours using a harmonic mean. They argue that the mixing length is constrained by the smaller of these two length scales.
To close this model we impose the mixing length to be equal to κ(h − d) at the top of the canopy (i.e. z = h, which is the bulk mixing length above an urban area in the standard roughness-length approach), as proposed by Coceal and Belcher (2004). This closure leads to the following formulation of the canyon mixing length l c : The displacement height d is determined by the empirical formulation proposed by Macdonald et al. (1998), which links the displacement height to the mean building height h and the building density λ p (often referred as the plan area index), which is defined as the ratio of the plan built area A plan to the total plan area A total : where α is an empirical parameter, whose chosen value is 4. Thus, the mixing length within the canopy is a function of morphological parameters of the canopy (h and λ p ). Finally, one can check that the model remains consistent with the extreme cases: -If the canopy is very sparse, then the density λ p tends toward 0, and so does the displacement height d. Thus, the mixing length tends towards the classical law of the wall (i.e. l m → κz), thereby reflecting the fact that the canopy does not impact the flow field.
-If the canopy is very dense, then λ p tends toward 1 and d ≈ h. Thus, l m tends toward l c and then the flow field is strongly influenced by buildings.

Wind profile
The Prandtl mixing model uses a logarithmic wind profile, which cannot be applied down to the ground in an urban canopy. Therefore, we instead used, within the urban area, an exponential profile, which is now widely used within vegetative canopies (Inoue, 1963;Petroff et al., 2008b). Numerous studies support the use of such a profile within the urban canopy (e.g. Macdonald, 2000;Masson, 2000). For example, measurements of median wind profiles within the urban canopy obtained during the Basel UrBan Boundary Layer Experiment (BUBBLE) are consistent with such an exponential wind profile within the urban canopy (Hamdi and Schayes, 2007). Assuming a mean flow above roof level, parallel to the canyon orientation, the exponential formulation is imposed all along the canyon. The exponential formulation can be deduced for a simple geometry (array of uniformly distributed drag elements), with simplifying hypotheses (mixing length and drag coefficient constant with height) as it was done for vegetative canopies: where β is an attenuation coefficient (Cionco, 1965) and u(h) the mean horizontal wind velocity at the building height h. Velocity profiles based on Eq. (16) are depicted in Fig. 4 for different values of β. One notes that, except for high values of β, the no-slip condition at the ground is not satisfied.
Based on studies by Arya (2001) and Rotach (1995), Masson (2000) computed the wind speed at half height for a narrow canyon (corresponding to the skimming flow). Subsequently, the following parametrisation of β was derived in this case: Hereafter, this expression will be assumed to apply for all canyon geometries. Another parametrisation of β is provided by Macdonald (2000), which is a linear relationship between the attenuation coefficient and the frontal building density λ f , defined as the ratio of the frontal built area A frontal to the total plan area A total : The sensitivity of the model to β is investigated in Sect. 5.4.4. The formulation, which was extracted from Masson (2000), was used in the following base simulations.
An integration over 360 • is performed to account for all street orientations. Only the wind component parallel to the canyon orientation is considered and thus a no mean wind condition inside the canyon is assumed when the flow is perpendicular to the canyon orientation: This formulation was computed for narrow urban canyons, i.e. for skimming flow conditions. Lemonsu et al. (2004) proposed to extend this formulation to all canyons. An adaptation of these formulations is used here. For wide canyons, in the case of the isolated roughness flow, the integration coefficient of the mean wind speed within the canyon is assumed to be equal to unity; subsequently, the formulation for wide canyons is the same as Eq. (16).
In the intermediate case, i.e. wake interference flow, the wind speed inside the canyon is computed as follows: We introduce for convenience sake the coefficient ζ , which depends on the canyon geometry, and we express the mean wind speed as follows: The no-slip condition requires that the wind velocity must be zero at the surface. Therefore, the exponential profile cannot apply near the surface and it must match with a different profile that tends to zero as z tends to zero. Experimental data suggest that, near the ground, the mean wind profile approaches a logarithmic profile (e.g. experimental data from Macdonald, 2000, Fig. 6).
The height z limit at which the change from the exponential wind profile to a logarithmic wind profile occurs is defined as the limit at which the mixing length in the urban canopy tends toward the law of the wall mixing length, i.e. l c κz limit l c + κz limit = (1 − )κz limit , where ∈ [0, 1] is a dimensionless parameter, which must be chosen as small as reasonably possible since too high value for will correspond to the assumption of a logarithmic profile for a large part of the urban canopy. The sensitivity of v d to the chosen value of is discussed in Sect. 5.4. The modified wind profile is depicted in Fig. 5.
The following formulation is assumed, according to the historical dry deposition velocity formulation (Gregory, 1945) where c(z ref ) is the concentration at the first vertical node of the air quality mesoscale model z ref (i.e. half the depth of the first model layer), v d is the dry deposition velocity seen from the atmosphere and F c atmosphere is the flux of pollutants removed from the atmosphere.
In order to compute the flux of pollutants removed from the atmosphere, the mass balance between the atmosphere and the surface can be written as follows, assuming there is no accumulation: where L is an area-averaged length of the street, defined by It should be noted that the canyon's width defines the exchange surface between the atmosphere and the canyon. These exchange surfaces are then defined at the top of the canopy between each region and the atmosphere.
Each F c canyon can be expressed by a mass balance in each region of the canyon: The different values of the dimensions of interest (fraction of street, wall and canyon which lie in the recirculation and the ventilation region) depend on the canyon geometry; they are summarised in Table 1. The term W wall refers to the height of walls. γ is defined, as the portion of the downwind wall, which lies in the recirculation region: We now describe the fluxes over each surface of the urban canyon.
Assuming Eq. (7), the mass flux can be written as In the case of turbulent mass transfer, the molecular diffusion term is negligible (aerodynamic resistance), whereas in the case of mass transfer in the quasi-laminar layer near the surface, the turbulent term is negligible (surface resistance).

Fluxes between the bulk atmosphere and the canyon
First, we assume here that the urban canopy is entirely contained within the first layer of the gridded air quality model. Second, we assume that the mass flux through the canyon is governed only by turbulent mass transfer. The flux from the bulk atmosphere (i.e. the atmosphere above the canyon) toward the canyon is chosen to occur from z ref to a reference height in the canyon region z canyon .
At this point, one must note that the well-known formulation of the dry deposition velocity depicted in Sect. 2 is based on the hypothesis of a constant vertical mass flux, which is not verified within the urban canopy, in particular the momentum flux formulation developed in this work is not consistent with this assumption (e.g. use of an exponential wind velocity profile). Nevertheless, in the absence of another available framework, we adapted this one-dimensional conceptual model of a vertical dry deposition flux to the twodimensional schematic representation of the urban canopy.
Accordingly, the flux is formulated as follows: with In the recirculation region, this integral is split into two parts, one above the canopy (z > h) and another one within the canopy (z < h). The continuity point is assumed to be the top of the canopy (z = h), as it was chosen for the improved formulation of the mixing length: Above the canopy (R top a, canyon ), the standard mixing length is used and the wind velocity is deduced from the classical logarithmic profile. The friction velocity u * is computed above the canopy with the Louis (1979) formula and parameters defining the canopy. The stability is then taken into account above the canopy. This integral can then be computed analytically. Within the canopy (R bottom a, canyon, r ), the improved mixing length is used (see Eqs. 13 and 14), and the wind velocity follows the exponential profile. This formulation leads to an indefinite integral (exponential integral Ei). It must be computed numerically.
In the ventilation region, the same resistance above the canopy is used (R top a, canyon ). Within the canopy, the mixing length, l m = κz, is used to reflect the weak influence of buildings on atmospheric turbulence in this part of the canyon.

Fluxes between the canyon and urban surfaces
For the sake of simplicity, in this section, the symbol means either street surface or wall surface, in each region (recirculation and ventilation). For the building walls and street surfaces, the flux can then be expressed similarly to the previous flux formulation where the concentration at the surface is taken to be zero. The turbulent mass flux occurs between z canyon and z 0, , which is the roughness length of the surface (building wall or street surface).
where R other, is either the surface resistance R s, in case of particles, or the sum of the quasi-laminar resistance and the surface resistance for gases, i.e. R b, + R c, .
In the recirculation region, the formulation of the aerodynamic resistance between the canyon and urban surfaces is expressed as follows: It is important to note that the roughness length z 0, represents the surface roughness and not the bulk roughness of the urban area. For the sake of simplicity, the aerodynamic resistance of the wall, is supposed to be similar to the aerodynamic resistance of the street, except for the local roughness length of the surface. A local friction velocity u * is also computed close to the surface. At this step, the atmospheric stability is assumed to be neutral.
The surface aerodynamic resistance in the ventilation region is written as follows: In the case when z limit is lower than z 0, , the logarithmic part of these equations is not taken into consideration and the lower bound of the integral is z 0, .

Flux between the bulk atmosphere and the building roofs
Dry deposition occurs also from the bulk atmosphere to the building roofs of the urban canyon area. The turbulence is assumed to be generated by the urban canopy above the roof. Subsequently, the following formulation applies: where R total, roof = R a,roof + R other,roof

Closure on the pollutant concentration within the canyon
The mass balance within the canyon (Eqs. 27 and 28) is used to close the flux equations and calculate the concentration c(z canyon ) needed in Eqs. (31) and (35) Note that, for the sake of simplicity, we have considered that no pollutant source is located within the urban canyon.

Overall dry deposition
The mass balance in Eq. (25)  atmosphere to an urban area:

Base simulation
It appears unfeasible to proceed to a quantitative comparison of the model proposed above to a set of measurements due to the paucity of dry deposition observations (see the discussion in Appendix B). However the model is applied to the Paris urban area, France, for the year 2011 and simulation results are compared to those obtained with a roughness-length model, the one described in Zhang et al. (2001). The meteorology is obtained from simulations conducted with the Weather Research and Forecasting Model (Skamarock et al., 2001). The surface resistances were computed following the model of Zhang et al. (2001), but the different local roughness lengths applied to walls and streets and the classical roughness-length approach apply to roofs lead to different surface resistances for these three types of surfaces. Calculations were performed here for particles with an aerodynamic diameter of 1 µm as an example. A single urban configuration is applied here over the whole domain for the sake of demonstration of the model; accordingly, a suburban configuration is assumed: -mean building height: h = 12 m; -mean roof width: W roof = 6.25 m (it is assumed that buildings are contiguous except for the side facing the street); -roughness length of walls: z 0,wall = 10 −4 m; -roughness length of streets: z 0,street = 10 −2 m. The dry deposition model presented above was implemented within the Polyphemus air quality modelling platform (Mallet et al., 2007). The roughness-length model based on Zhang et al. (2001) was already available in the Polyphemus platform. The meteorological fields are interpolated from the Weather Research and Forecasting (WRF) discretization grid to the Polyphemus grid. After this preprocessing, meteorological data are provided with a horizontal resolution of 0.04 • × 0.027 • every hour. For land use coverage, the Global Land Cover 2000 database (http://bioval.jrc. ec.europa.eu/products/glc2000/glc2000.php) is used and the 23 original categories are aggregated to match the land use categories defined by Zhang et al. (2001). Outside urban areas, the roughness-length model based on Zhang et al. (2001) was used. Figure 8 shows the mean dry deposition velocity computed with the parametrisation presented in this work (λ p = 0.4). These results are consistent with the range of measurements reported in the literature (see Appendix B). Figure 9 represents the mean wind speed at the reference height z ref averaged over the year 2011. The dry deposition velocity is strongly influenced by the mean wind speed, inasmuch as the aerodynamic resistance depends on it. Greater deposition velocities occur in areas with greater wind speeds. Figure 10 shows the annual average over the year 2011 of the hourly relative difference between the dry deposition velocity computed with the model presented above (λ p = 0.4), v canyon , and computed with the roughness-length model, Geosci. Model Dev., 8, 893-910, 2015 www.geosci-model-dev.net/8/893/2015/ The differences are computed for each hour then averaged over the year 2011. The mean over all the fully urban grid cells (100 % of urban coverage) of the annual average of the hourly relative differences is about 45 % with a mean standard deviation (SD) of 18 % (not shown). This mean difference reaches 82 % for λ p = 0.6 with a SD of 26 % (not shown). These relatively low SDs are explained by the fact that the two models use similar approaches, based on wind velocity profiles. The different vertical wind profiles, mixing length and surface areas used in the two models explain the difference. In Fig. 11, the time series of this relative difference during a winter period (from January to March 2011) is presented for different building densities for one chose urban grid cell. The sensitivity of the deposition velocity to the value of the building density is discussed in Sect. 5.3.

Total flux over urban surfaces
A major difference between the standard roughness-length model and the model developed here, is the ability of the latter to distinguish surfaces within the urban canopy. Figure 12 depicts the normalised dry deposition rates over walls (black crosses and black line), roads (red crosses and red line) and roofs (blue line), which have been calculated for the Paris suburbs (λ p = 0.4) in November 2011 (a uniform pollutant concentration of 1 µg m −3 is used to normalise the deposition rate). Figure 13 depicts the dry deposition fluxes on each surface and region for the same period. These fluxes are also normalised with a unit atmospheric concentration of 1 µg m −3 .
The major fraction of dry deposition fluxes occurs on the roofs. The resistance to deposition is strongly influenced by Relative difference (%) ∆v d (roughness, λ p =0.4) ∆v d (roughness, λ p =0.6) ∆v d (roughness, λ p =0.8) Figure 11. Time series of dry deposition velocity relative difference between the urban canopy and roughness-length models from January 2011 to March 2011 in one chosen urban grid cell for three values of the building density. the distance to the surface; thus, the deposition flux on roofs is larger than on any other surface.
In this configuration, even with a building density, λ p = 0.4, the ventilation region is very narrow, and its contribution is close to zero, even if fluxes on surfaces in this region are significant (see Fig. 13). It explains the reason why the deposition rate is close to zero in this region (strictly zero for the street, because there is no portion of the street that lies in the ventilation region). Figure 13 also shows that the modelled deposition on building walls is slightly lower than on streets in the same region (there is no sedimentation on walls). In the present parametrisation, the modelled deposition fluxes in the ventilation region are slightly larger than in the recirculation region. This difference can be explained by the rather strong shear layer created in the recirculation region, which implies that this part of the canyon is nearly isolated from the bulk atmosphere and explains that the deposition fluxes are lower in this region (see Fig. 13).

Influence of building density
The impact of the building density (λ p ) on the dry deposition velocity was investigated. In Fig. ,  λ p = 0.8, which is a rather theoretical density.
The dry deposition velocity computed with the roughnesslength model is slightly lower than the dry deposition velocity computed with the present parametrisation for low to medium mean wind speed. At high wind speed, the dry deposition velocity computed with the roughness-length model crosses over the one computed with the urban canopy model for a very low building density (λ p = 0.2). The dry deposition velocity computed with the roughness-length model is nearly linear with the wind speed, whereas the one computed with the urban canopy model is not.
The difference between the roughness-length model and the urban canopy model can be partly explained by the fact that the surface available for deposition is greater in the lat- ter. However, when the building density is very low, additional deposition surfaces are not large enough to compensate the resistance of the last few metres, which are not taken into account in the roughness-length model (i.e. from d to the ground).
Concerning the urban canopy model at higher building densities, one might expect that, as the turbulence increases, the deposition rate should grow with building density. However, once a threshold is exceeded (λ p ≈ 0.6), the dry deposition velocity decreases with the building density in the present parametrisation, because the strong shear layer generated by buildings nearly suppresses interactions between the canyon and the bulk atmosphere (i.e. R a,canyon,r increases strongly). Such results can only be obtained with an urban canopy model that provides some differentiation among different flow regimes within urban canyons. These results are consistent with measurements obtained in Chicago and South Haven by Yi et al. (2001). They found that the dry deposition velocity (overall dry deposition velocities for various pollutants) was higher in Chicago (moderate λ p ) than in South Haven (low λ p ).

Other sensitivity tests
In this section, the sensitivity of the model results to the following key parameters is investigated: the coefficient α of the displacement height formulation, the characteristic recirculation length W r , the threshold z limit and the attenuation coefficient β of the exponential wind profile. Macdonald et al. (1998) chose to set the coefficient α in Eq. (15) to 4. This value was obtained from experiments conducted over arrays of cubes. We have tested our model, using the following canyon characteristic lengths: As it can be seen in Fig. 15, the dry deposition velocity is not very sensitive to the values of α. For λ p = 0.8 and α ∈ [2, 6], the dry deposition velocity varies by less than 7 %. It varies by less than 6 % for λ p = 0.6 and by 4 % for λ p = 0.4. Therefore, the canopy model is not very sensitive to the value of α and the default value of 4 seems appropriate.

Characteristic recirculation length W r
The canyon characteristic recirculation lengths are defined empirically as 2 to 3.5 times the building height (see Sect. 3.1) Therefore, for a building height of 12 m, we conducted simulations with W r varying from 24 to 42 m. For λ p = 0.6 and λ p = 0.8, the dry deposition velocity is not very sensitive to the value of W r (not shown). It can be explained by the fact that, in these densely built configurations, the ventilation region is very narrow or nonexistent. For a rather low building density (λ p = 0.4), which includes a large ventilation region, the dry deposition varies by only 8 % (not shown).
Therefore, this parameter has little influence on the canopy model. It can affect the distribution of pollutants within the canopy because it defines the boundary between the recirculation region and the ventilation region, but it has little effect on the amount of pollutants removed from the atmosphere.

Threshold z limit
The threshold z limit defines the height at which the wind profile within the urban canyon switches from a logarithmic pro- file near the surface to an exponential profile. The sensitivity of the dry deposition velocity to the threshold z limit (via the value of ) is investigated. For all building densities, the variation does not exceed 9 % for ∈ [0.1, 0.2] (not shown). The dry deposition velocity is not very sensitive to the value of z limit . A value of = 0.2 was chosen.

Attenuation coefficient β of the exponential wind profile
The sensitivity of the dry deposition velocity to the attenuation coefficient β is illustrated in Fig. 16. This parameter strongly influences the dry deposition velocity. According to Cionco (1972), β should be between 0.5 and 3 for a wide range of vegetative canopies. In this range, the dry deposition velocity varies by a factor of about 2.
Several formulations are available to define β (see Sect. 3.2.2). In order to decide which formulation should be chosen, it is important to note that the geometry chosen in this work (λ p = λ f ) does not seem to be compatible with the ones studied by Macdonald (2000). In fact, the frontal building density in his work is considered to be lower than 0.35. Above this density, the formulation could not be applied. Since the building density in this work varies between 0.2 and 0.8, MacDonald's formulation was not considered here. Regarding Masson's formulation, it is within the range recommended by Cionco (1972). Moreover, it has been computed from measurements in a real urban area (Toulouse, France), and then confirmed with another experiment (e.g. Lemonsu et al., 2004, during the ESCOMPTE campaign). Therefore, Masson's formulation was used here.

Conclusions
The standard roughness-length model is appropriate if one is interested in the removal of airborne pollutants from the atmosphere. However, if one wants to follow the spatial distribution of pollutant deposition within urban areas, the roughness-length model is not suitable because it fails to differentiate among the different types of surfaces (roofs, walls, streets, etc.). For example, the experimental results of Roupsard et al. (2013) suggest that dry deposition velocities can vary by a factor of 24 between two surface types in urban areas. Consequently, there is a need to be able to model dry deposition in urban areas with some spatial resolution.
We have presented an urban canopy model for dry deposition that takes into account the atmospheric flow regimes depending on urban morphology and resolves three types of surfaces (roofs, walls and streets). Therefore, this approach provides three-dimensional spatially distributed dry deposition fields within the urban canopy, which cannot be obtained from the roughness-length model. The model was shown not to be very sensitive to key parameters related to the atmospheric flow within the urban canyon (except for the attenuation coefficient β). The building density affects the dry deposition velocity. For a suburban area, the urban canopy model led to greater dry deposition than the roughness-length model. For sparsely built areas, both modelling approaches gave similar results at low wind speeds but diverge at high wind speeds due to their different vertical wind profile formulations. For very densely built areas, the formation of a shear layer prevents dry deposition within the urban canyon and there is an optimal building density that maximises dry deposition in the present model. This work has shown that fluxes of pollutants may vary by a factor of 4 among different surfaces and regions in a given urban area. Further work could address a finer characterisation of surface materials in terms of roughness to better estimate dry deposition fluxes according to surface types, thereby leading to even greater variability in dry deposition fluxes.
Further work could also address finer representations of micrometeorology within the urban canopy (e.g. improved wind profile). Using a meteorological model with fine vertical resolution (or a multi-layer model such as the one proposed by Martilli et al., 2002) within the urban canopy could provide valuable information on atmospheric turbulence. An adaptation of a multi-layer canopy model could also be used to refine the aerodynamic resistance formulation. Furthermore, applications to actual urban configurations should be conducted.
Above the urban canopy, the thermal stability has been taken into account through the use of the Louis (1979) parametrisation. Within the urban canopy, a neutral condition was assumed. Because of the urban heat island, the layer within the urban canopy could be either neutral or unstable. Further analysis should be conducted to investigate the influence of unstable conditions on the dry deposition flux.
The contribution of sources of some air pollutants (e.g. from road traffic) within the urban canopy may need to be taken into account.
Finally, there is a need for measurements of dry deposition velocities in urban areas, which could be used to evaluate and improve dry deposition models. However, such measurements are very difficult to implement and new method developments (e.g. use of isotopes in laboratory settings) may be required to obtain experimental databases suitable for model performance evaluation. For particles, the sedimentation velocity must be added to the mass transfer by diffusion. Therefore, the particle mass flux is expressed as follows: where v s is the gravitational settling velocity. The turbulent mass flux through the canyon surface for particles can be formulated as follows, under the same hypothesis as for gases: 1 − e (−v s R a, canyon ) .
Likewise, for a street surface, the mass flux formulation for particles can be expressed as follows: For a wall surface, since gravitational settling is not relevant for mass transfer to vertical walls, F c wall can be expressed with the same equation as that for gases (Eq. 35).

Appendix B: Overview of dry deposition observations
There is a wide range of existing methods to measure dry deposition velocities and it is of interest to discuss briefly the advantages and limitations of such dry deposition velocity measurements (Wesely and Hicks, 2000). There are two main quantification methods of deposition: direct and indirect measurements. Direct measurements use surrogate surfaces that mimic the actual surface, and are used to quantify dry deposition via analysis of the amount of material deposited on the surface. However, although the use of surrogate surfaces is convenient for the collection and analysis of material, it raises the question of representativeness of such surfaces compared to actual surfaces. Moreover, it has been shown that both the surface geometry and characteristics have a large impact on deposition (Sakata and Marumoto, 2004).
Indirect measurements are typically based on micrometeorological approaches where the dry deposition flux is calculated by measuring both the atmospheric concentration and the vertical velocity. There is a wide range of techniques to measure fluxes (eddy correlation, eddy accumulation, gradient method, etc.). These techniques provide a flux measurement that is representative of a large homogeneous area. Even the interpretation of these measurements remains questionable (see Baldocchi et al., 2000, for instance). Thus, such techniques cannot provide detailed information on dry deposition fluxes in complex settings such as urban areas. Furthermore, there are very few experimental data available on dry deposition over urban areas. The scarce field campaigns generally occur over long sampling period of time and detailed meteorological data over these periods are not available.
Such studies are extremely difficult to conduct because of the heterogeneity of the environment, large spatial and temporal variations of meteorological conditions and the challenges associated with the measurements of dry deposition fluxes. One can cite dry deposition fluxes measured around Lake Michigan Paode et al., 1998;Shahin et al., 2000;Yi et al., 2001, for example) during the Atmospheric Exchange Over Lakes and Oceans Study (AE-OLOS). Yi et al. (2001) measured overall dry deposition velocities in Chicago, which vary from 2.1 cm s −1 (fine particle fraction of Cu and Zn) to 23 cm s −1 (fine particle fraction of Al and Mn). In another study, Sakata and Marumoto (2004) measured dry deposition on the roof of a building at Komae City, Japan. They measured dry deposition velocities in a range from 0.73 cm s −1 (for Zn) to 4.6 cm s −1 (for Mn). Moreover, these measurements show a high standard deviation, which makes interpretation difficult. Sofuoglu et al. (1998) reported a factor of 5 between measured and modelled (using a roughness dry deposition velocity) fluxes in Chicago.
Clearly, the large uncertainties associated with dry deposition flux measurements make their use for model performance evaluation difficult, as dry deposition can vary by more than 1 order of magnitude depending on surface type and meteorological conditions (e.g. see Roupsard, 2013 2 , for an exhaustive review of dry deposition velocity measurements over urban surfaces).

Code availability
The code supporting this study can be found in the Supplement of the manuscript.
The Supplement related to this article is available online at doi:10.5194/gmd-8-893-2015-supplement.