The efficient urban canopy dependency parametrization (SURY) v1.0 for atmospheric modelling: description and application with the COSMO-CLM model for a Belgian summer

This paper presents the Semi-empirical URban canopY parametrization (SURY) v1.0, which bridges the gap between bulk urban land-surface schemes and explicitcanyon schemes. Based on detailed observational studies, modelling experiments and available parameter inventories, it offers a robust translation of urban canopy parameters – containing the three-dimensional information – into bulk parameters. As a result, it brings canopy-dependent urban physics to existing bulk urban land-surface schemes of atmospheric models. At the same time, SURY preserves a low computational cost of bulk schemes for efficient numerical weather prediction and climate modelling at the convectionpermitting scales. It offers versatility and consistency for employing both urban canopy parameters from bottom-up inventories and bulk parameters from top-down estimates. SURY is tested for Belgium at 2.8 km resolution with the COSMO-CLM model (v5.0_clm6) that is extended with the bulk urban land-surface scheme TERRA_URB (v2.0). The model reproduces very well the urban heat islands observed from in situ urban-climate observations, satellite imagery and tower observations, which is in contrast to the original COSMO-CLM model without an urban land-surface scheme. As an application of SURY, the sensitivity of atmospheric modelling with the COSMO-CLM model is addressed for the urban canopy parameter ranges from the local climate zones of http://WUDAPT.org. City-scale effects are found in modelling the land-surface temperatures, air temperatures and associated urban heat islands. Recommendations are formulated for more precise urban atmospheric modelling at the convection-permitting scales. It is concluded that urban canopy parametrizations including SURY, combined with the deployment of the WUDAPT urban database platform and advancements in atmospheric modelling systems, are essential.


Introduction
Cities over the world are expanding (Seto et al., 2012) and an increasing share of the population tends to live in the cities (United Nations, 2014).The associated changes to the landscape and anthropogenic heating led to excess temperatures in cities compared to their natural surroundings, which is known as the urban heat island (UHI) effect.The UHI causes a higher exposure to heat stress leading to excess mortality (Laaidi et al., 2011;Gabriel and Endlicher, 2011), damage to infrastructure, higher energy usage (e.g.indoor active cooling) intensifying outdoor urban heat and greenhouse gas emissions and pressure on socio-economic activities.In view of global climate change with more extreme heat waves Published by Copernicus Publications on behalf of the European Geosciences Union.ahead, vulnerabilities and impacts are increasing for urban centres of all sizes, economic conditions and site characteristics across the world (Revi et al., 2014).
During the past 3 decades, a vast amount of urban landsurface schemes have been developed.They enable the convection-permitting atmospheric models to resolve the heterogeneity of cities with applications for heat stress assessment and the development of urban climate adaptation and mitigation strategies (Prein et al., 2015).Even though their purpose of representing urban physics in land-surface schemes of atmospheric models is the same, intercomparison studies (Karsisto et al., 2016;Trusilova et al., 2016;Best and Grimmond, 2015;Grimmond et al., 2011) demonstrate that they differ in terms of modelling strategy, complexity, input parameters and applicability: on the one hand, the bulk schemes (e.g.Meng, 2015;De Ridder et al., 2015;Wouters et al., 2015;Pleim et al., 2014;Grossman-Clarke et al., 2005;Fortuniak et al., 2004) take into account the overall radiative and thermal properties, vegetation sparseness, the surface roughness, the water-storage capacity and the anthropogenic activity in the urban canopy with a set of bulk parameters.These model parameters are estimated from model sensitivity experiments (Wouters et al., 2015;De Ridder et al., 2012;Demuzere et al., 2008;Fortuniak, 2005), observational campaigns (Masson et al., 2008;Rotach et al., 2005;Offerle et al., 2005) and inventories (Flanner, 2009;Pigeon et al., 2007).The bulk schemes are suitable for capturing the general urban climate characteristics at the convection-permitting scales.They include the interplay between the excess conversion of incoming radiation into sensible heat, the heat accumulation, reduced wind speeds and the additional anthropogenic heating, which results in the urban heat island effect (Phelan et al., 2015;Zhao et al., 2014).However, they do not explicitly resolve the complex processes depending on the local characteristics and the three-dimensional structure of the urban canopy, which further modulate the urban climate.These processes include shadowing and multiple scattering of radiation, heterogeneous surface-atmospheric interaction in terms of turbulent momentum, heat and moisture transport and the innerbuilding energy budget.On the other hand, the more complex explicit-canyon schemes (Demuzere et al., 2014;Allegrini et al., 2014;Schubert et al., 2012;Hénon et al., 2012;Oleson et al., 2008;Fortuniak, 2007;Kanda et al., 2005;Martilli et al., 2002;Masson, 2000;Crawley et al., 2000) explicitly capture one or more of these complex physical processes.As a result, they allow for representing the detailed micro-scale features, such as heterogeneous temperatures of the facets and canyon wind gusts.Therefore, their added value compared to the less complex bulk schemes lies in their ability to acquire the more detailed information about urban climate risk and vulnerability at the micro-scales.Moreover, Best and Grimmond (2015) have shown that an improved skill in offline model experiments is found when increasing the number of explicitly resolved urban processes.Yet, the ap-plicability of these explicit-canyon schemes for convectionpermitting atmospheric modelling is generally hindered by either the lack of detailed urban canopy information, computational cost and their model complexity, and these issues are explained in more detail hereafter.
First, the complex schemes are hindered by the lack of urban canopy parameters, which include information about three-dimensional urban morphology and material properties.Detailed parameter inventories are available for specific urban sites.On the one hand, they are suitable for extensive offline evaluations of the urban land-surface schemes (e.g.Loridan and Grimmond, 2012;Grimmond et al., 2011).On the other hand, such detailed canopy inventories form a basis for spatially varying data sets with a world-wide coverage, such as Faroux et al. (2013), Loveland et al. (2010), Jackson et al. (2010) and Bartholomé and Belward (2005).These data sets are intended to account for the urban canopy parameter variability at the global scales.Furthermore, more detailed spatially varying data sets exist, but they only include specific parameters over specific cities.An example here is the use of morphological data from CityGML (Gröger and Plümer, 2012) in urban climate modelling (Schubert and Grossman-Clarke, 2014).Despite those initiatives, expansion and updates of urban canopy information remain challenging (Seto et al., 2011) at the convection-permitting scales, not to mention the prognoses from future land-use change scenarios (Prein et al., 2015).Consequently, the spatial detail, accuracy, coverage and variety in those urban canopy parameters data sets are limited, which is also clear from the substantial differences between the data sets (Schneider et al., 2009).However, the first urban model intercomparison project demonstrates that such parameter information is important for improved modelling performance in existing urban land-surface schemes (Best and Grimmond, 2015).For those complex schemes, missing data could deteriorate the model performance (Grimmond et al., 2011).In order to tackle these issues, the WUDAPT (World Urban Database Access Portal Tools) initiative has been developed recently, aiming for a coherent and detailed urban canopy parameter data set with a world-wide coverage on a 100 m resolution supported by a peer network of researchers.It allows the urban-climate research for taking into account the urban canopy variation in residential, commercial and industrial areas.However, the database is under development and consequently its coverage is still limited.
Second, explicit-canyon schemes are computationally demanding compared to bulk schemes.Complex schemes could lead up to 15 % of total computational cost of an atmospheric model (e.g.Trusilova et al., 2016).In contrast, bulk schemes allow very fast downscaling of ensemble climate projections (e.g.Lauwaet et al., 2015).Finally, the complexity of the explicit-canyon schemes sometimes makes a consistent incorporation into the host atmospheric model and its maintenance very challenging.This paper presents the Semi-empirical URban canopYparametrization (SURY) (Sect.2.1) which accomplishes the following: 1.It offers a translation of urban canopy parameters containing the three-dimensional information into bulk parameters.The translation is based on detailed observational studies, modelling experiments and available parameter inventories.
2. It brings canopy-dependent urban physics -which used to be reserved for the explicit-canyon schemes -to existing urban bulk urban land-surface schemes (e.g.Meng, 2015;De Ridder et al., 2015;Wouters et al., 2015;Pleim et al., 2014;Grossman-Clarke et al., 2005;Fortuniak et al., 2004).Therefore, it allows them to take into account the city's heterogeneity in the urban canopy parameters.
3. It preserves the low computation cost and low complexity of the bulk schemes.
4. It provides versatility in employing either urban canopy parameters from bottom-up inventories or bulk parameters from top-down estimates.
As a result, SURY bridges the gap between bulk schemes and explicit-canyon schemes by providing robust canopydependent urban physics and parameter versatility for convection-permitting atmospheric modelling at a low computational cost.
SURY is implemented in the COSMO(-CLM) model (Sect.2.2).The latter is extended with the bulk urban landsurface scheme TERRA_URB v2.0, which allows for taking bulk parameters from SURY into account.The model system is set up over Belgium during a mid-summer period in 2012 at 2.8 km resolution (Sect.2.3).It is compared with observational data from in situ urban-climate observations of air temperature and land-surface temperatures from satellite imagery.As shown by Wouters et al. (2013) with an idealized boundary-layer advection model, the nocturnal boundary-layer urban heat island (BLUHI) intensity and its vertical extent, reaching their maximum during the night, are affected by the vertical temperature profile in the lowest few hundreds of metres above the ground.Therefore, an evaluation is also performed against the nocturnal boundary-layer temperatures from tower observations.The meteorological measurements are described in Sect.2.4.The COSMO(-CLM) model without an urban land-surface scheme has already been evaluated extensively for numerical weather prediction at the convection-permitting scale (Baldauf et al., 2011) and regional climate modelling applications (Brisson et al., 2016b, a;Bucchignani et al., 2016;Prein et al., 2015;Thiery et al., 2015;Vanden Broucke et al., 2015;Fosser et al., 2014;Ban et al., 2014;Feldhoff et al., 2014;Dosio et al., 2014;Kotlarski et al., 2014;Lange et al., 2014;Van Weverberg et al., 2014;Panitz et al., 2013;Berg et al., 2012;Keuler et al., 2012, and references herein).Therefore, the evaluation (Sect.3.1) focusses on urban-climate modelling, particularly the thermal contrast between cities and the natural surroundings.It is investigated whether the model can reproduce, the canopy-layer UHI (CLUHI; estimated from the urban-rural differences in screen-level temperature), the surface UHI (SUHI; urban-rural differences in the land-surface temperatures) and the nocturnal boundary-layer urban heat island (BLUHI; urban-rural differences in the vertical temperature profiles).As an application of SURY, the effect of urban canopy parameter changes on regional climate modelling with the COSMO(-CLM) model is quantified with an online sensitivity experiment (Sect.3.2).The sensitivity allows for addressing the possible effect of the parameters' uncertainty and variability on urban-climate modelling, hence allowing for setting priorities in acquiring them.Finally, a discussion and concluding remarks are given in Sect. 4.

Semi-empirical urban canopy parametrization
In this section, the SURY is described.The translation of urban canopy parameters into urban bulk parameters is based on the urban physical processes with regard to the ground heat transfer, the surface-radiation exchanges and the surface-layer turbulent transport for momentum, heat and moisture: the bulk thermal parameter values take into account enhanced ground heat transport and storage due to the increased contact surface with the atmosphere (Fortuniak et al., 2004) expressed by the surface area index (SAI) in Sect.2.1.1.Furthermore, the radiative bulk parameter values take into account the albedo reduction factor resulting from the radiative trapping by the urban canopy in Sect.2.1.2.Finally, the enhanced surface drag on the wind by the buildings in the urban canopy take into account the building height in Sect.2.1.3.As a result, SURY introduces an efficient dependency of bulk urban land-surface schemes on the urban canopy parameters.Throughout the subsections below, the robustness of SURY is verified by comparing bulk parameters from top-down estimates with those translated from bottom-up urban canopy parameter inventories.In addition, the default urban canopy parameters are obtained from the recommended values for the medium urban category in Loridan and Grimmond (2012;see their Table 4, stage 5b), for which the respective bulk parameters are determined.An overview of the urban canopy parameters (SURY input) and the bulk parameters (SURY output) is given in Table 1.

Ground heat transport
A new methodology is developed for translating the urban canopy parameters into bulk (effective) thermal parameters for heat capacity and heat conductivity.The latter is taken into account in the slab representation, which considers the www.geosci-model-dev.net/9/3027/2016/Geosci.Model Dev., 9, 3027-3054, 2016 one-dimensional heat equation of a vertical column commonly used in existing land-surface schemes.In the methodology, the buildings and pavements are considered as massive impermeable structures stacked on the natural soil.It takes into account the three-dimensional surface curvature of the urban canopy, which results in a larger contact surface with the atmosphere than a slab surface enhancing the ground heat flux.As denoted by Fortuniak et al. (2004), this results in generally larger urban bulk thermal admittances, µ bulk = λ bulk C v,bulk (De Ridder et al., 2012;Demuzere et al., 2008), than the corresponding urban material values (Loridan and Grimmond, 2012;Jackson et al., 2010).The bulk values of the heat capacity and heat conductivity at the surface level are obtained by multiplying the corresponding urban canopy values with the SAI.Therefore, SAI is the ratio between the land-surface area that envelops the urban canopy and the plan area.In this way, the ground heat transport is integrated over the land-surface area.Such a multiplication is substantiated in more detail as follows.The vertical surface heat transport through the uppermost layer is calculated by Fourier's law: where λ s (unit: W m −1 K −1 ) is the surface heat conductivity and ∂T ∂z is the vertical temperature gradient.At the same time, the tendency of the vertical profile is described by the heat equation (assuming that the upper-layer heat conductivity is constant with depth): where C v,s is the surface heat capacity (unit: J m −3 K −1 ).We now consider the case that the surface exposed to the atmosphere is enlarged with the SAI factor, while conserving the original vertical temperature profile of the uppermost layer.On the one hand, more heat goes through the larger surface, hence the total heat flux through the surface is multiplied with the same SAI factor.In order to get such an enhanced heat transport in the slab representation, λ s needs to be multiplied with the SAI factor as well.On the other hand, the tendency of the vertical temperature profile of the new uppermost layer remains the same as the original (this will be the case if the atmospheric forcing would not change, which is true for an infinitesimal time period after enlarging the surface).For the latter, a heat equation for the slab representation is required that is equivalent to the original equation (i.e.Eq. 2), hence λ s /C v,s ratio needs to be unchanged.In conclusion, the multiplication of the ground heat flux by a factor SAI and the conservation of the original ground temperature tendency are required for the slab representation.These requirements can be attained with the original set of onedimensional equations and using bulk parameters for which both λ s and C v,s are multiplied with the same factor SAI.
In the case of an idealized urban canopy with parallel urban canyons, straight roads and flat roofs, SAI can be ob-tained from geometrical considerations: where h w c is the canyon height-to-width ratio and R the roof fraction.The first right-hand term of this equation represents the surface area index of the street canyon.In turn, it is subdivided in 1 × (1 − R) which is the surface area of the street, and 2 h w c ×(1−R) which is the surface area of the two walls in the street canyon.By adding the roof fraction R, one obtains the total surface area index of the idealized urban canopy.
A methodology is presented that simultaneously takes into account the bulk surface thermal properties of the urban canopy (derived above) and those of the natural soil below the urban canopy.Therefore, vertical profiles of the bulk thermal parameters C v,bulk (z) and λ bulk (z) are calculated, with z the vertical coordinate in the bulk model (ground depth).As explained above, the bulk surface heat capacity is obtained by multiplying the surface heat capacity C v,s with the surface area index SAI: Below the surface, the urban substrate layer with a thickness equal to the building height h is considered for representing the thermal mass of the urban canopy in thermal contact with the natural soil below.The bulk heat capacity in this layer considers a vertical linear gradient between the surface value (C v,bulk,s ) and the value of the natural soil below (C v,soil ): Below the urban substrate layer, the bulk heat capacity is equal to C v,soil : An analogous formulation is considered for the vertical profile of the bulk heat conductivity λ bulk (z): where λ soil is the heat conductivity of the natural soil and λ bulk,s is the bulk surface heat conductivity: The default urban canopy parameters (see also Table 1) are set equal to the recommended values for the medium urban category in Loridan and Grimmond (2012); see their Table 4 (stage 5b).Herein, the height of the buildings h, the canyon height-to-width ratio h w c and the roof fraction R are equal to 15 m, 1.5 and 0.667, respectively.According to Eq. ( 3), the values for R and h w c led to an SAI of 2.0.The default values for the surface heat conductivity C v,s (1.25 × 10 6 J m −3 K −1 ) and the surface heat capacity λ s (0.767 W m −1 K −1 ) are the respective weighted averages from the values for roof ).Therefore, the weighted averages are calculated according to the surface fractions of roofs, walls and roads in the urban canopy: where the first terms represent the contributions from the street canyon and the second terms those from the roofs.The values for C v,bulk,s and λ bulk,s are obtained from Eqs. ( 4) and ( 9) and they yield 2.5×10 6 J m −3 K −1 and 1.53 W m −1 K −1 , respectively.
The bulk surface thermal admittance is expressed as Given the values for C v,bulk,s and λ bulk,s above, one obtains µ bulk,s = 1.96 × 10 3 J m −2 K −1 s −1/2 .It lies within the range for the thermal admittance of the compact and open climate zones in Table 4 of Stewart and Oke (2012) and also within the uncertainty range obtained by De Ridder et al. (2012).Although this is not a formal validation, these correspondences give confidence to the default parameter values of C v,s and λ s above and to the enhanced effective surface heat capacity and heat conductivity expressed by Eqs. ( 4) and (9).It needs to be noted that the presented methodology above assumes a homogeneous surface temperature of the urban canopy.This is also case for the next section with regard to the surface radiation properties.Consequentially, the scheme does not explicitly represent the temperature variety among the different elements in the urban canopy resulting from shadowing and the heterogeneous thermal and radiative properties.Therefore, urban-physical processes resulting from such variety are not explicitly resolved.This choice was made for providing consistency with the bulk urban landsurface schemes employing bulk parameters.

Surface radiation
In this section, the methodology for deriving the bulk (or effective) albedo α bulk and emissivity bulk from urban canopy parameters is addressed.The bulk values refer to the portions of reflected incoming short-wave radiation and emitted infra-red radiation by the urban canopy layer to the upper atmosphere, respectively.It also accounts for the modulation of the bulk value according to the increased-albedo effect of snow.The bulk albedo reduction factor of the urban canopy ψ bulk is derived from the canyon height-to-width ratio where α is the surface albedo and f snow is the snow-covered fraction.ψ bulk h w c , R is calculated by where ψ c ( h w c ) is the canyon albedo reduction factor.Therefore, the total albedo reduction factor is calculated from the albedo reduction of the roof weighted with the roof fraction R and that of the street canyon weighted with the canyon fraction (1 − R).As stated before, flat roofs are considered, hence the roof fraction R does not lead to a albedo reduction.In contrast, multiple reflections take place for the street canyon for which the canyon albedo reduction factor ψ c is taken into account.
Instead of implementing a computationally demanding explicit-canyon radiation scheme, an approximation for ψ c is proposed to the numerical estimation from Fortuniak (2007).The latter applies an exact solution of the multiple-reflection problem allowing to subdivide the different facets in an urban canyon.The exact solution results in a high accuracy for low solar heights when the lower canyon parts are shaded.It could reproduce the effective-albedo observations from a scale model (Aida, 1982) and from a real canyon very well.The numerical estimation shows that the albedo reduction is most sensitive to the h w c ratio, hence the following approximation is proposed: This closely matches the numerical estimation with a maximal error of ±7 % for the highest excursion of the sun during summer solstice at the mid-latitude (55 • ), a canyon parallel to the solar azimuth and an albedo of 0.4 (Fortuniak, 2007; see their Fig.11).With regard to other sun heights, canyon directions relative to the solar azimuth and h w c ratios between 0 and 2 (Fortuniak, 2007; see their Fig.8), the proposed ψ c formulation has a maximal error of 45 %.It should be noted that the approximation is fitted to the numerical estimation for a perfect urban canyon.Hence, the approximation neglects additional albedo changes due to bending roofs and varying albedos for the different facets.
Optionally, a distinction is made between the albedo of roofs, roads and walls as follows: with and where is the averaged albedo of the roads and walls in the urban canyon.The bulk infra-red emissivity bulk takes into account the same bulk albedo reduction factor ψ bulk as follows: where is the emissivity and snow is the snow emissivity.
The robustness of the Eqs.( 15), ( 16) and ( 17) is verified for a dense urban area in Toulouse centre: the averaged albedo for walls, roofs and roads are 0.25, 0.15 and 0.08, respectively, whereas the roof fraction and h w c ratio are 0.59 and 1.4, respectively; see Pigeon et al. (2008).This yields a snow-free bulk albedo α bulk for the urban canopy of 0.125.This is close to the bulk value of 0.13 for summer, which is estimated from the averaged ratio between the upward and the downward radiation measured by a mast tower during the CAPITOUL campaign (Masson et al., 2008).
The default urban canopy parameters (see also Table 1) for albedo of roofs (0.10), walls (0.10) and roads (0.15) are set equal to the recommended values from Loridan and Grimmond (2012); see their Table IV (stage 5b).Together with the values for h w c = 1.5 and R = 0.667, Eq. ( 16) yields a (snow-free) bulk albedo of α bulk = 0.081.At the same time, the bulk albedo reduction factor for the urban canopy yields ψ bulk = 0.80 (see Eq. 14).With respect to the more simple formulation in Eq. ( 13), a surface albedo of α = 0.10 is obtained and used by default.Analogously, the default values for the snow-free bulk emissivity of bulk = 0.89 and the surface emissivity of = 0.86 are obtained.

Surface-layer turbulent transport
Following Sarkar and De Ridder (2010), the aerodynamic roughness lengths for the urban canopy is calculated as follows: with h as the building height.The thermal roughness length z 0H is obtained with a parametrization of the inverse Stanton number (as in De Ridder, 2006;Demuzere et al., 2008): Geosci.Model Dev., 9, 3027-3054, 2016 www.geosci-model-dev.net/9/3027/2016/ with k as the von Kàrmàn constant.For the urban canopy, a bluff-body thermal roughness length parametrization from Brutsaert (1982) is introduced using parameter values from Kanda et al. (2007): where Re * = u * z 0 /ν is the roughness Reynolds number, u * is the friction velocity and ν = 1.461×10 −5 m 2 s −1 the kinematic viscosity of air.
As before, the default value for h equal to 15 m (see also Table 1) corresponds to the recommended value in Loridan and Grimmond (2012); see their Table 4 (stage 5b).It yields z 0 = 1.125 m and kB −1 = 13.2 (in case that u * = 0.25 m s −1 ).

The COSMO(-CLM) model
The COSMO model (Steppeler et al., 2003) is a full 3-D atmospheric numerical model designed for operational and research applications in limited-area weather prediction at high resolution.The model has been developed by the German Weather Service (DWD) and is further improved and maintained by the Consortium for Small-Scale Modelling (COSMO).Members of this consortium include several meteorological services from inside and outside Europe.
The COSMO model has a compressible non-hydrostatic core for atmospheric dynamics and includes parametrizations for radiative transfer, cloud microphysics, subgridscale turbulence and convection.It also includes parametrizations for the ground heat and water transport and the landatmosphere interactions, as described in more detail in the next paragraph.The regional climate model COSMO-CLM (COSMO model in CLimate Mode) is based on the COSMO model and includes modifications allowing the application on timescales up to centuries (Böhm et al., 2006;Rockel et al., 2008).These modifications comprise the introduction of an annual cycle to vegetation parameters like the vegetation cover and the leaf area index as well as an externally prescribed, time-dependent CO 2 concentration in the atmosphere.The COSMO-CLM model is further developed by a vast amount of researchers of the CLM community inside and outside of Europe (http://www.clm-community.eu).It is used extensively for long-term regional climate studies (Klutse et al., 2015;Endris et al., 2015;Vanden Broucke et al., 2015;Cavicchia et al., 2014;Davin et al., 2014;Schubert and Grossman-Clarke, 2013) and for downscaling global-climate realizations (Akkermans et al., 2014;Dosio and Panitz, 2015;Lejeune et al., 2015).
The ground heat and water transport and the representation of vegetation and snow cover are resolved by the soil-vegetation-atmosphere transfer (SVAT) module TERRA_ML (Schulz et al., 2016;Doms et al., 2011;Grasselt, 2008).The COSMO(-CLM) model implements the next-generation TKE-based surface-layer transfer scheme (Doms et al., 2011;Buzzi, 2008).The surface layer, which refers to the layer of air between the earth surface and the lowest model level, is divided into a laminar-turbulent sublayer, the roughness layer and a constant-flux (or Prandtl) layer.The surface layer scheme is also intimately related to the TKE-based closures of the COSMO(-CLM) model; see Sect. 3.3 and 3.4 of Doms et al. (2011).As a result, the surface layer does not need the empirical Monin-Obukhov stability functions (as in Paulson, 1970;Guo and Zhang, 2007) for which the Obukhov stability parameter needs to be determined from an iterative procedure or a non-iterative approximation (e.g. Louis, 1979;Wouters et al., 2012;Li et al., 2014, and references herein).It rather generates these functions by the use of the dimensionless coefficients of the turbulence closure (Mellor and Yamada, 1982); see Sect.4.2 of Doms et al. (2011).Land-surface parameters including the soil type, vegetation and orography are specified with the external parameter tool (EXTPAR); see also Smiatek et al. (2008).These are processed from global land cover data sets, such as those for orography and coastlines from Global 30 Arc-Second Elevation (GTOPO30; see https://lta.cr.usgs.gov/GTOPO30)or ASTER (Three Arc-Second Elevation data set), soil data from the Digital Soil Map of the World (DSMW; see http://www.fao.org/geonetwork/srv/en/metadata.show?id=14116) and land-use data from Global Land Cover 2000 (GLC2000; see Bartholomé and Belward, 2005), GLOBCOVER (Loveland et al., 2010) or ECO-CLIMAP (Faroux et al., 2013).The vegetation parameters, which include vegetation cover fraction, LAI and rooting depth, are specified with annual minimum and maximum values depending on the land use (Doms et al., 2011; see their Tables 14.3, 14.4 and 14.5).For the extra-tropical Northern Hemisphere, a growing and resting period is calculated according to the latitude.The roughness length parameter (over land) depends on both land use and the subgrid-scale orography.The obtained vegetation and soil data are assigned to the natural soil fraction.The methodology for the reduction in vegetation abundance in urban environments compared to the rural surroundings is according to the underlying landuse data set.
In the original COSMO(-CLM) model, cities are represented by natural land surfaces with an increased surface roughness length and a reduced vegetation cover.However, in this representation, urban areas are still treated as waterpermeable soil with aerodynamic, radiative and thermal parameters similar to the surrounding natural land.Therefore, this basic representation could not reliably capture the urban physics and associated urban-climatic effects including urban heat islands.In order to tackle this issue, the bulk scheme TERRA_URB (Wouters et al., 2015) has been introduced for providing an intrinsic representation of urban physics in the COSMO(-CLM) model.The modified ground heat and moisture transport and land-atmospheric exchanges for momentum, heat and moisture found over urban areas are adopted by including modifications to the soil-vegetation module TERRA_ML (Grasselt, 2008;Schulz et al., 2016), the land-atmosphere interactions (Doms et al., 2011) and the input data (EXTPAR; Smiatek et al., 2008).The initial release (version 1) features a non-iterative calculation of surface-layer stability functions accounting for the roughness sublayer (Wouters et al., 2012) and an impervious waterstorage parametrization based on a probability density function of water puddles (Wouters et al., 2015).TERRA_URB v1.0 has been evaluated in offline mode in Wouters et al. (2015) for intensive observation campaigns in Basel (Rotach et al., 2005) and Toulouse (Masson et al., 2008).It has also been employed for acquiring heat-stress scenarios of future climate change and urban land-use change scenarios in Belgium (De Ridder et al., 2015), adopted by the Climate Report of the Flemish Environmental Agency (Brouwers et al., 2015).During the Online Urban Model Intercomparison Project (Trusilova et al., 2016), TERRA_URB v.1 has been compared to other urban land-surface parametrizations TEB (Trusilova et al., 2013;Masson, 2000) and DCEP-BEP (Schubert et al., 2012;Martilli et al., 2002) coupled to same COSMO-CLM model.The next version (v2.0) of TERRA_URB introduces several advancements compared to its previous version.The main advancements are the implementation of SURY v1.0 and the application of the TKE-based surface-layer turbulent transfer module of the COSMO(-CLM) model (Doms et al., 2011;Buzzi, 2008).A full description of TERRA_URB v2.0 can be found in Appendix A.

Model setup
The COSMO(-CLM) model that implements SURY in its urban land-surface module TERRA_URB v2.0 is set up for the reference simulation over Belgium; see Fig. 1.The simulation is performed during a summer period from 1 July to 20 August 2012 for which the first 3 weeks are considered as spin-up.The model parameter setup is based on the COSMO-CLM model configurations of Brisson et al. (2016b, a) and Prein et al. (2015), employing convection-resolving climate simulations.The domain covers an area of 175 × 175 grid cells centred over Brussels with a horizontal grid spacing of 2.8 km resolution.A total of 40 vertical layers are used with the lowest domain level at 10 m above the ground.For the lateral boundaries, the model takes 3-hourly analysis from the operational model of the European Centre for Mediumrange Weather Forecasting (ECMWF) at a spatial resolution of 0.125 • in latitude and longitude.The reference simulation above, referred to as "REF", is compared with a simulation ("STD") that uses the original COSMO-CLM model without urban land-surface parametrization.The surface parameters for TERRA_URB v2.0 are the spatially detailed impervious surface area (ISA) data set from the environmental agency (Maucha et al., 2010) and a data set of anthropogenic heat emission (AHE) inventory modified from Flanner (2009); see Appendix A6.As indicated in the introduction, existing global databases for the urban canopy parameters are not suitable for the intended scale of cities and neighbourhoods.Even though the first steps have been made by mapping its local climate zones (Verdonck et al., 2016), more detailed inventories and databases, specifically from http://WUDAPT.org, do not cover yet the Belgian region.As an intermediate solution, the default urban canopy parameters from Sect.2.1 are employed (see upper panel of Table 1), which contemplates the combination of available parameter inventories, modelling and observational studies and SURY's theoretical framework.
In addition to the REF and STD setup described above, a range of online-coupled sensitivity experiments are performed.Starting from REF for each sensitivity simulation, the parameters for the urban canopy are changed according to the minimum (low scenarios; L) and maximum values (high scenarios; H) of the urban canopy parameter ranges derived from the local climate zones of compact low-rise and Except for the AHE, L and H correspond to the minimum and maximum values of the urban canopy parameter ranges for the local climate zones of compact low-rise and mid-rise defined in Stewart and Oke (2012).For the GL scenario, the AHE is set to 0 W m −2 .For the GH scenario, AHE multiplied by 2 compared to the default setup for which the data set and methodology of Flanner (2009, denoted  mid-rise in Stewart and Oke (2012).The parameters are, respectively, the surface albedo (α), surface heat conductivity (λ s ), surface heat capacity (C v,s ), canyon height-to-width ratio h w c , building height (h) and roof fraction (R).The parameter range for α is obtained from the range in bulk albedo α bulk in Stewart and Oke (2012), while keeping the other (morphological) parameters at their default values.Therefore, Eq. ( 13) is reversed to calculate α from α bulk .In the same way, the urban canopy parameter ranges for the surface heat capacity and heat conductivity are derived from the range of the bulk thermal admittance in Stewart and Oke (2012).In addition, the AHE is set to 0 W m −2 in the L scenario and it is multiplied by 2 in the H scenario.These two scenarios are in agreement with the given uncertainty range in Stewart and Oke (2012) -indicating a range between 0 and 75 W m −2 -while preserving daily and annual variability in the model and spatial variability from the input data.As a result, 14 additional simulations are performed for which the L scenarios are, respectively, AL, BL, CL, DL, EL, FL and GL, and the H scenarios are AH, BH, CH, DH, EH, FH and GH.An overview of the simulations can be found in Table 2.

In situ measurements
The modelled screen-level air temperature and the associated CLUHI is evaluated against in situ measurements for an urban and a rural site in Antwerp.These were performed using platinum resistance thermometers supplied by Campbell Scientific.The sensors were mounted in an actively ventilated radiation shield (Young 43503) to reduce heating effects by radiation loading on the sensor and stagnant air inside the shield.The sensor plus radiation shield setups were deployed side by side during a test phase during 1 week at the end of June 2012.The root mean square difference on the 15 min temperature averages during this week was found lower than 0.04 • C. The rural station was located on the premises of an organic farm enterprise 10 km to the south-east of Antwerp.The station is situated in pastureland with grass kept short by sheep.The nearest buildings are low rise and about 200 m away.The urban station was located on the premises of the Royal Lyceum in Antwerp, a secondary school.The sensor was mounted in the urban canopy layer on top of the roof of a 3 m high small building in the centre of a fairly large playground.The distance to the nearest adjacent wall was about 15 m.Though the location on top the roof is probably not ideal in terms of micro-scale effects, comparisons with similar measurements in the vicinity of the lyceum and an analysis of the urban-rural differences suggest that the measurements at the lyceum can be considered representative for the neighbourhood.The position of the in situ measurements are indicated in Fig. 1.

MODIS satellite imagery
As demonstrated recently by Hu et al. (2014) and Tomlinson et al. (2012), satellite information can be used for evaluating the models' excess surface heating in urban areas compared to the natural surrounding areas for the entire study domain.The model results are evaluated against land-surface temperatures (LSTs) derived from the MODIS sensor on board of the Terra and Aqua satellites, providing four overpasses a day with a spatial resolution of approximately 1 km.Hu et al. (2014) compared three methods to evaluate MODIS LST to HRLDAS surface temperatures and the optimal technique is employed for this study.For each Terra and Aqua overpass, the nearest model time is selected resulting in images at 09:00 and 11:00 UTC and 20:00 and 00:00 UTC.The first two are combined as "day", while the latter two are referred to as "night".Cloudy pixels are removed in MOD11A1 and MYD11A1 version 5 (Wan, 2008), while potential remaining cloudy pixels are removed by only using MODIS LST data within 1.5 times the interquartile range (Hu et al., 2014;Monaghan et al., 2014).In addition, cloudy pixels from the COSMO-CLM model are removed in both data sets.Finally,

Tower observations
In order to evaluate the nocturnal boundary-layer temperature and BLUHI in the model, observations from two meteorological towers within the province of Antwerp (Belgium) are used.The first tower of the Flemish Environmental Agency (VMM) is 160 m high and located on an industrial site in Zwijndrecht (geographical coordinates: 51 • 14 37.9 N, 4 • 20 3.1 E).Measurements for temperature, humidity, wind speed, pressure and precipitation are performed at 8, 24, 48, 80, 114 and 153 m.The temperature measurements are obtained with multi-stage solid state thermistor (Met One 062).The second tower of the Belgian Nuclear Research Centre (SCK•CEN) is 120 m high and located in Mol (geographical coordinates: 51 • 13 4 N, 5 • 5 24 E).Measurements have been done in the framework of determining atmospheric stability according to the turbulence scheme of SCK•CEN (Bultynck and Malet, 1972).It is placed in a rural area on flat terrain with pine trees in the immediate vicinity.Measurements for wind direction, wind speed and temperatures are performed at 8, 24, 48, 69, 78 and 114 m above ground level.Therefore, the temperature mea-surements are obtained with a Cu-constantan thermocouples.The temperature measurements of both towers are placed within a ventilator-driven aspirated radiation shield to protect against direct and diffuse solar heating.They are continuously recorded on a 1 min basis in order to make and store 10 min averages.For the model evaluation, the available temperature measurements at the heights 8, 48 and 114 m of both towers have been used.The positions of the tower observations are indicated in Fig. 1.

Evaluation
The averaged LST and screen-level temperatures for the dayand night-time from the REF simulation during the summer evaluation period are shown in Fig. 2. It is found that the SUHI reaching its maximum during the day is typically larger (approximately 4 K for the urban centres) than the CLUHI reaching its maximum during the night (approximately 2 K).The night-time CLUHI is of comparable magnitude as the night-time SUHI.The CLUHI reaches a its minimum at noon.Therefore, the urban heat islands occur at the scale of the cities and magnitude increases with city size.These findings are consistent with existing literature for urban-climate modelling and observational studies (e.g.Phelan et al., 2015).Therefore, excess conversion of solar radiation into sensible heat is taking place in cities at daytime.The excess heat is partly stored into the urban canopy which leads to day-time excess land-surface temperatures.In turn, the prolonged excess heat release from the urban canopy substantially affects the shallow nocturnal boundary layer resulting in augmented screen-level temperatures at night-time (Wouters et al., 2013;Bohnenstengel et al., 2011).Modelled temperatures and associated urban heat islands of the REF simulation are evaluated in more detail below against screen-level temperatures from the in situ measure-ments (Sect.3.1.1),land-surface temperatures (LST) from satellite imagery (Sect.3.1.2)and nocturnal boundary-layer temperature from tower observations (Sect.3.1.3).Herein, a comparison is also made with the STD simulation that excludes the urban parametrization.

Two-metre air temperatures
The model evaluation of the screen-level air temperature and the CLUHI for Antwerp is displayed in Fig. 3   For the urban site, reference model REF agrees well with the observed temperatures with a mean bias (MB) of +0.39 K and a Pearson correlation coefficient (R 2 ) of 0.96.In contrast, the standard model STD results in a general temperature underestimation with a MB of −0.95 K and an R 2 of 0.95.For the rural site, the temperature time series daily variability are also well reproduced for both REF and STD with an MB 0.56 and 0.14 K and an R 2 of both 0.95, respectively.However, the model shows an underestimation of the diurnal cycle, particularly an overestimation of the nocturnal temperatures.The overestimation for the rural site is larger for REF than for STD.It is postulated that the latter stems from the ex- The CLUHI intensity (1.74 K), calculated from the observed difference between the urban station and the rural station, is well reproduced by REF (1.56 K) with a very good R 2 of 0.80 and a bias of −0.18 K.In contrast, the CLUHI could not be reproduced by the STD model; with an averaged UHI value of 0.64 K, STD has a much larger model bias −1.10 K and a lower R 2 of 0.54.Again for the REF simulation, an underestimation of the UHI occurs during the night and a slight overestimation in the late afternoon.The modelled standard deviation of 1.25 K is lower than the observed 1.56 K, which originates from the underestimation of the nocturnal UHI.
For the NU class, a negative bias of −1.53 K is found for the day-time LST samples and a positive bias of 4.36 K for the night-time samples for both REF and STD.The too-low diurnal cycle is consistent with the results for the screen-level air temperatures (see Sect. 3.1.1).The model results for REF confirm the observed gradual increase of LST with the increasing urban density.At the same time, the day-and nighttime SUHI increases as well with the urban density in both REF and the observations.For the day, the modelled SUHI by REF has a negative median bias of −1.52 and −1.14 K for the MU and DU class, respectively, whereas it is −0.6 K for the LU class.For the night, the averaged modelled SUHI by REF matches very well the observations.Therefore, the biases in SUHI for the LU, MU and DU classes are 0.03, 0.26 and 1.09 K, respectively.This general positive nocturnal bias is in contrast to the negative bias for the screen-level air temperatures (see Sect. 3.1.1).In contrast to REF, STD is not able to capture the SUHI.

Nocturnal boundary layer
The model results are evaluated with observations of the nocturnal boundary-layer described in Sect.2.4.3.Besides absolute temperature profiles, we also focus on the performance of the nocturnal BLUHI.The latter is calculated from the difference between temperature profiles in an industrial and natural area.The results are shown in Fig. 5 and Table 6.As before, sensitivity results are included, which will be addressed in Sect.3.2.Because of their height, it needs to be noted that profile differences between the two mast towers also originate from other sources that are not related to the urbanization.In particular for certain wind directions, substantial influence could arise from the Scheldt River near the tower of Zwijndrecht.
The REF simulation is capable of capturing the temperature profiles and the variability of both towers.For the industrial site, the model profile matches well that of the observations with a positive bias of 0.41 K.A higher positive temperature bias is found for the rural site (0.73 K), which generally stems from the temperature overestimation near the ground.Therefore, the increasing temperature with height ( stable profile) in the model output is overestimated.The overestimated temperatures near the ground are consistent with the positive nocturnal bias in the screen-level measurements shown in Sect.able to reproduce the contrast between the urban and rural site, particularly the more stable boundary layer in the rural site compared to the industrial site.As a result, the observed boundary-layer UHI effect (1.02 K), calculated as the difference between the temperature profiles of the industrial and the natural site, is well reproduced by REF (0.70 K) with a negative bias of −0.32 K and a correlation of 0.58.The positive bias for STD (0.31 K) for the rural tower is slightly lower than that of REF.As similar as in Sect.3.1.1,it is postulated that the latter stems from the excess urban heat introduced by the urban land-surface scheme in REF, which is advected from the urban areas towards the rural areas.A much larger temperature shift between STD and REF occurs for the industrial tower (ZWN).The boundary-layer temperatures in STD are underestimated by −0.86 K, which is in contrast to REF for which it is overestimated.As the tendency towards more neutral/unstable temperature slopes in the industrial site is absent in STD, the boundary-layer UHI could not be reproduced with a negative bias of −1.17 K and a correlation 0.29.

Urban parameter sensitivity
The model sensitivity experiment in response to the different urban input parameters is performed by means of modelbias analysis.Therefore, the model output and performance statistics in terms of the observed temperature quantities (described in Sect.2.4) for the different experiments (AL-GH) described in Sect.2.3 are compared with the reference simulation (REF) and the simulation without urban parametrization (STD).The results of the sensitivity experiments can be found in Figs. 3, 4 and 5 and in Tables 3, 4, 5 and 6 (Taylor diagrams are not shown).In general, the different model sensitivities are restricted to the urban areas and therefore they emerge as changes in the UHIs.The model sensitivity is much smaller than the model changes when including urban parametrization (REF minus STD).This is the case for screen-level temperatures, LST and nocturnal boundary- With respect to the screen-level temperature and intensity for Antwerp, the sensitivity range is generally larger during the evening and during the night when the CLUHI remains close to its maximum.The largest sensitivity range for the diurnal cycle is found for the parameters AHE and (to a smaller extent) the thermal parameters, which are the surface heat conductivity and heat capacity.A medium sensitivity range to the diurnal cycle is found for the building height, canyon height-to-width ratio and the roof fraction.The overall lowest model sensitivity is found from the albedo, yet at day-time, when overall sensitivity changes are the lowest, it is higher than the sensitivity from the thermal parameters.For the daily averaged CLUHI, the magnitudes of the different sensitivities are comparable.A more detailed description of the different parameter sensitivities is found below: -An increased AHE obviously leads to an increase in screen-level urban temperatures and heat island intensity.Even though the AHE is higher during the day, the sensitivity range from AHE is the largest during the night because of the confinement in the nocturnal boundary layer (Wouters et al., 2013;Bohnenstengel et al., 2011).
-A diurnal change in sign of the sensitivity is found for the thermal parameters: lower thermal parameter values led to a temperature increase during the day followed by a decrease during the night.This is due to the fact that a smaller portion of excess heat is needed to heat up the urban canopy during the day, hence leading to higher day-time temperatures.At the same time, less excess heat buffering in the urban canopy occurs which lowers the temperatures at night-time.The sign change results in a relatively low sensitivity of the daily-mean temperature, even though the day-and night-time sensitivities are relatively high.
-An overall decrease in day-and night-time CLUHI appears from the decrease in building height.The consequential lower roughness length (Eq.19) results in higher wind speeds.As a result, the urban excess heat is less accumulated in the urban centres lowering the air temperatures there.Smaller buildings also yield a lower effective heat capacity and heat conductivity of the soil layers below (see Eqs. 5-8).In the same way as in the previous point concerning the thermal parameters, such a change counteracts (yet not fully) the abovementioned day-time temperature decrease and enhances the night-time temperature decrease.
-Both the increase in roof fraction and surface albedo lead to an overall decrease in nocturnal and day-time UHI, because they both increase effective albedo according to Eq. ( 13).It is remarkable that the sensitivity www.geosci-model-dev.net/9/3027/2016/Geosci.Model Dev., 9, 3027-3054, 2016 to the roof fraction is larger than that for the albedo at night-time, whereas this becomes reversed at day-time.This means that another process different from the effective albedo effect is also important; in fact, a larger roof fraction also implies a reduction in effective heat capacity and heat conductivity (see Eqs. 3, 4 and 9).Again, in the same way as in the previous point concerning the thermal parameters, this counteracts the decrease in day-time UHI, while it enhances the decrease in night-time UHI.
-A lower canyon height-to-width ratio h w c leads to an overall decrease in screen-level urban temperatures and UHI for the day and for the night.Two mechanisms are at play here: on the one hand, a decrease in short-wave radiative trapping in the canyon (see Eqs. 13 and 14) leads to a higher effective albedo, hence a reduced conversion of solar radiation into heat.On the other hand, a lower h w c ratio decreases the heat transfer below the surface due to the lower contact surface (see Eq. 3) implying a lower effective heat capacity and conductivity (see Eqs. 4 and 9).Combined, they lead to an additional decrease in nocturnal urban temperature and UHI, while counteracting the day-time decrease resulting from the higher effective albedo mentioned above.
The sensitivity results for LST sensitivity are consistent with the sensitivity to the screen-level temperatures.For example, lower surface heat conductivity or heat capacity leads to higher urban LSTs and SUHI at day-time and lower urban LSTs and SUHI at night-time.However, the hierarchy of sensitivity to LST and SUHI is different from that of the screen-level temperature: the largest sensitivity range for the urban LST and SUHI now relates to the thermal parameters and not to the change in AHE, which now has a medium sensitivity range.This is explained by the fact that AHE is considered as a heat source to the first atmospheric model layer (see Appendix A4), hence only indirectly influencing the LST.Sensitivity of LST to the other urban canopy parameters are similar to those of the screen-level temperatures: the h w c ratio and roof fraction yield a medium sensitivity range, and the building height yields a low sensitivity range to LST.The sensitivity to the albedo is low for the night-time LST and high for the day-time LST.The sensitivity results for the temperature profiles show a similar sensitivity as for the screen-level temperatures.Therefore, the model sensitivity propagates throughout the nocturnal boundary layer for at least the first 110 m above ground level, hence it modulates the BLUHI.
Analysis of the model biases and Taylor plots for each of the temperature quantities (not shown) demonstrate that the model performance change for each of the parameter sensitivity runs is ambiguous and depends on the temperature quantity considered.On the one hand, parameter changes leading to a better performance in CLUHI and BLUHI sometimes lead to worse performance in terms of absolute temperatures.In particular, a decrease in the value of the thermal parameters leads to a lower positive model bias and better model performance in terms of absolute screen-level temperatures, but a larger negative bias and worse model performance in terms of CLUHI.On the other hand, the perfor- mance change between the SUHI and CLUHI may lead to a divergent behaviour as well: an overall improvement in the CLUHI (lower positive day-time bias; lower negative nighttime bias) by increasing the value of the thermal parameters leads to an overall deterioration of the SUHI (higher day-time negative bias; higher night-time positive bias).Finally, the day-and night-time performance changes are also divergent.For example, increasing the albedo leads to an overall deteriorated performance in day-time LST and SUHI, but leads to an improved night-time performance.

Discussion and conclusions
SURY is developed for efficiently representing canopydependent urban physics in atmospheric models.The methodology bridges the gap between bulk and explicitcanyon schemes in atmospheric models.The urban canopy parameters are translated into bulk parameters.The latter can be easily taken into account in bulk urban land-surface schemes in existing atmospheric models.The urban canopy parameters as input for SURY include the canyon heightto-width ratio, the building height, the roof fraction, the short-wave albedo, thermal emissivity and the heat conductivity and heat capacity.The output parameters of SURY include bulk albedo, bulk emissivity, aerodynamic and thermal roughness length and vertical profiles of the bulk heat conductivity and heat capacity.The methodology delivers theoretical and empirically verified robustness that is based on detailed urban observational and modelling experiments.Additional model robustness has been provided by comparing existing bulk parameters from top-down estimates with those translated from bottom-up urban canopy parameter inventories.
The outline of SURY implemented in an atmospheric model system is given in Fig. 6.SURY, in combination with any existing bulk urban land-surface scheme, provides several benefits over other methodologies for urban canopy parametrization in atmospheric modelling at the convectionpermitting scales.They are related to applicability, computational cost and model consistency: -SURY allows to employ detailed bottom-up data sets representing the variation in residential, commercial and industrial areas, as soon as they are available from WUDAPT or any other newly available data set.Still, one can fall back on the bulk parameters when the detailed information is missing.Such a parameter versatility enhances the overall applicability of the bulk landsurface schemes.It also enables the analysis of the propagation of uncertainty in the input parameters by comparing simulations using bulk parameter estimates with those using urban canopy parameter data sets.
-SURY v1.0 enables the comparison and consolidation between bottom-up urban canopy parameter data sets (translating them to bulk parameters) and top-down bulk parameter data sets.This provides a way to verify consistency of the parameter data sets.It should be noted that such a parameter consolidation has been done in Sect.2.1 on the basis of existing parameter inventories leading to the parameter list in Table 1.
-SURY combined with a bulk urban land-surface parametrization is beneficial for efficient convectionpermitting atmospheric modelling applications.In particular, the bulk albedo approximation avoids explicit numerical computation of the complex canyon radiation trapping, hence largely reduces the computational demand compared to other explicit-canyon radiation modwww.geosci-model-dev.net/9/3027/2016/Geosci.Model Dev., 9, 3027-3054, 2016 els.In case of the COSMO(-CLM) model, one only measures a computational overhead of 7 % over the original COSMO(-CLM) model without urban landsurface parametrization.The majority of the computational overhead stems from the implementation of the tile approach in TERRA_URB and the overhead from SURY itself is negligible.
-As SURY translates urban canopy parameters into bulk parameters, it can be easily applied for intrinsic implementation of urban physical processes and their dependency on urban canopy parameters in existing atmospheric models already employing a bulk urban landsurface scheme.An intrinsic implementation instead of an external urban land-surface parametrization preserves consistency between the representation of the different physical processes and features in urban environments on the one hand and those of natural environments on the other hand.Therefore, it is assured that the contrasts in model response between urban areas and natural areas only stems from the differential landsurface parameter values.In this way, artificial discrepancies that arise from possible different formulations between an external land-surface module resolving the urban areas and an internal land-surface module resolving the natural areas are avoided.Furthermore, it allows for the transparent application of intrinsic features of the host atmospheric model in the urban-physics modelling, such as the representation of vegetation (shading) in the urban canopy.In the case of the COSMO(-CLM) model, these features include the multi-layer snow representation and the TKE-based surface-layer turbulent transfer scheme.
SURY is evaluated in online mode with the COSMO(-CLM) model.To this end, the COSMO(-CLM) model is extended with the bulk urban land-surface scheme TERRA_URB v2.0, which allows for taking urban bulk parameters from SURY v1.0 into account.An online model evaluation over the Belgian urban extent during summer 2012 demonstrates that the model is capable of capturing urban climate characteristics.The model captures well the daily and diurnal variability of the UHIs in terms of landsurface temperatures, screen-level temperatures and nocturnal boundary-layer temperatures.Therefore, most of the negative temperature biases occurring for the urban areas in the original COSMO(-CLM) model without urban parametrization are alleviated.
Although TERRA_URB v2.0 implementing SURY v1.0 offers a general model improvement with regard to urban temperatures and urban heat islands, it could not alleviate other systematic errors in the COSMO(-CLM) model.Such systematic errors include an overall underestimation of the diurnal cycle of absolute temperatures, particularly an overestimation of the nocturnal temperatures.This is in agreement with previous evaluations (regarding summer) of the COSMO-CLM model (Brisson et al., 2016b;Trusilova et al., 2016;Brisson et al., 2016a;Lange et al., 2014;Keuler et al., 2012).This is especially the case for the rural absolute temperatures for which the urban effect is small and can be related to the underestimated stability of the nocturnal temperature profile.One should keep in mind that the systematic errors could be propagated, amplified or compensated by the urban-atmospheric interactions, hence deteriorating urbanclimate modelling (assessment).Particularly, an overall underestimation is found for the diurnal cycle of the UHI in agreement with Trusilova et al. (2016).Therefore, an underestimation of the nocturnal UHI is found, which may originate from the overestimation of nocturnal temperatures and underestimation of stable temperature profiles found in rural areas.This can be expected from the fact that a reduced stability of the nocturnal boundary-layer temperature profile of the city leads to a reduction of the canopy-layer UHI, as demonstrated with an idealized boundary-layer advection model in Wouters et al. (2013).
The model sensitivity is investigated for which SURY takes urban canopy parameter ranges of the local climate zones of compact low-rise and compact mid-rise in Stewart and Oke (2012).The model response and performance change to urban canopy parameter changes are generally restricted to the urban areas.They are also smaller than the effect of the introduction of the urban land-surface scheme in the atmospheric model.With regard to the screen-level temperatures, the nocturnal boundary-layer temperature and associated CLUHI and BLUHI, the largest model sensitivity is found for the AHE and the urban thermal parameters, a medium sensitivity for the building hight and h w c ratio and the lowest sensitivity for roof fraction.With regard to the dayand night-time LST and SUHI, the largest model sensitivity is found for the urban thermal parameters, a low sensitivity for the AHE, h w c and roof fraction and a very low sensitivity for the building height.For both SUHI and CLUHI, the sensitivity to the albedo is relatively high yet slightly lower than for the thermal parameters during day-time, whereas it is slightly lower than AHE, h w c and roof fraction during nighttime.
The following recommendations are made with regard to future urban-climate research: -To date, limitations in either accuracy, detail, variety or coverage exist in urban canopy parameter inventories.(Stewart and Oke, 2012).Therefore, Bechtel et al. (2015) provide an easy way to perform such classification (level 0 information), which has recently been applied for the three largest cities in Belgium (Verdonck et al., 2016).The combination between such LCZ mapping and local urban canopy parameter sampling and inventories forms a basis for providing the spatial detail in the urban canopy parameters (See et al., 2015).Such urban canopy parameters can be directly used as input for SURY (see upper panel of Table 1) or any other urban canopy parametrization.In addition, the level 1 and 2 data collection, providing information on the form and function within each LCZ class, can deliver additional parameter detail.
-As SURY allows for verifying consistency between urban canopy parameters and bulk parameters, its application is recommended in the WUDAPT framework.In turn, this allows urban-climate research for more precise climate assessment.Moreover, it will lead to more consistency in the comparison of bulk schemes with explicit-canyon schemes in future urban model intercomparison projects.
-The model performance largely depends on the temperature quantity considered.In particular, parameter settings leading to better UHI sometimes lead to worse absolute temperatures and vice versa, which is also the case for day-time temperatures vs. night-time temperatures, and land-surface temperatures vs. air temperatures.This ambiguity demonstrates that a multi-variable model evaluation is a requirement for improving and comparing urban-climate modelling strategies.
-Some of the aforementioned model errors exceed the model sensitivity range with regard to the urban canopy parameter uncertainty.This demonstrates that the majority of the model uncertainty may not be related to urban canopy parameter uncertainty, but may result from deficiencies in the land-surface module and other aspects of the coupled atmospheric model.Such systematic errors may also cause the ambiguity in model performance and its sensitivity to the urban parameters.In addition to the advancements of urban canopy parametrizations and the expansion of more detailed parameter databases, ongoing improvements in atmospheric modelling systems are also essential for more precise urban-climate modelling assessment.Regarding the representation of urban heat islands, focus should be given to those components that improve the representation of the surface fields and boundary-layer characteristics, such as the stability of the nocturnal temperature profiles.
-In view of the advantages listed above, consideration of employing SURY is recommended for numerical weather prediction and long-term regional climate modelling applications at the convection-permitting scales.
Besides the presented implementation in the COSMO-CLM model, this can also be achieved by applying SURY in other existing urban bulk land-surface schemes for intrinsic implementation in atmospheric models.Herein, it should be kept in mind that SURY's semi-explicit nature implies some limitations with respect to complex urban physical processes in an urban environment.Particularly, the heterogeneity of the urban canopy induces micro-scale dynamic and physical processes in the urban canopy, such as shadowing, buoyancy flows and the complex inner-building heat and moisture transfer.As SURY does not resolve the full heterogeneity, the propagating micro-scale features -such as heterogeneous temperatures patterns at the street level and canyon wind gusts -are not represented.
In case there is need for a more detailed urban-climate impact assessment that explicitly accounts for such features, one requires more detailed urban-physics modelling with higher complexity.They include the more explicit radiation schemes allowing heterogeneous temperatures between sun-lit and shaded facets (e.g.Schubert et al., 2012), computational fluid dynamics (e.g.Allegrini et al., 2014) and inner-building energy models (e.g.Bueno et al., 2012;Crawley et al., 2000).It should be noted that these methodologies require much finer meshes, larger computational resources and very detailed information on urban design.
It is concluded that urban canopy parametrizations including SURY, combined with the deployment of the WUDAPT urban database platform and advancements in atmospheric modelling systems, are essential.They allow for more prewww.geosci-model-dev.net/9/3027/2016/Geosci.Model Dev., 9, 3027-3054, 2016 Appendix A: TERRA_URB v2.0 This appendix describes TERRA_URB v2.0, the bulk urban land-surface scheme of the COSMO(-CLM) model.Modifications to the surface-layer turbulent transport scheme of the COSMO(-CLM) model, the modified surface evaporation and transpiration, the anthropogenic heat emission, the tile approach and the additional surface input parameters are successively given in the subsections below.Bulk parameters are calculated from the urban canopy parameters with the SURY described in Sect.2.1; see also Table 1.

A1 Surface-layer turbulent transport
In contrast to its previous version employing the similaritybased turbulent transfer scheme of Monin and Obukhov (1954), TERRA_URB v2.0 makes use of the TKE-based surface-layer turbulent transfer scheme of the COSMO(-CLM) model.Herein, the exchanges for momentum and heat within the surface layer between the surface and the lowest atmospheric model layer are determined as follows: where ρ and c p are the density and the specific heat of the air, u A and θ A are the absolute wind speed and potential temperature at the lowest atmospheric model layer, θ S is the potential temperature of the surface, u * is the friction velocity and θ * is the turbulent temperature scale.Stability-dependent transfer-layer resistances (i.e.resistances to the surface-layer turbulent transfer between the surface and the lowest atmospheric model layer) for momentum (r M SA ) and heat (r H SA ) are calculated with the TKE-based surface-layer transfer scheme of the COSMO(-CLM) model (Doms et al., 2011).For the urban canopy, aerodynamic roughness length is obtained from Eq. ( 19).For the natural land cover, z 0 is adopted from standard input parameters of the COSMO(-CLM) model; see Sect.2.2.The thermal roughness length z 0H is obtained with a parametrization of the inverse Stanton number (as in De Ridder, 2006;Demuzere et al., 2008): ) with k as the von Kàrmàn constant.For the natural areas, the inverse Stanton number is diagnosed from the TKE-based surface-layer transfer scheme as follows: where r H S0 is the roughness-layer resistance for heat (i.e. the resistance to the surface-layer turbulent heat transfer between the surface and z 0 ), e = 0.5(u 2 + v 2 + w 2 ) the turbulent kinetic energy and S H (z 0 ) the stability function for the Prandtl layer.For the urban canopy, the bluff-body thermal roughness length is used from Eq. ( 21).Inserting Eq. ( 21) into Eq.(A4) leads to a modified r H S0 .In turn, this enters the surface-layer heat transfer equation (Eq.A2) as follows: where r H 0A is the free atmospheric resistance (i.e. the resistance between the surface-layer turbulent transfer z 0 and the lowest atmospheric model layer).

A2 Evaporation and transpiration
In contrast to its previous version only considering the water on the impervious surfaces, TERRA_URB v2.0 is also capable of explicitly representing the water-permeable surfaces (bare soil), vegetation and snow in the urban canopy.The evapotranspiration is obtained in a similar way as the natural land in the standard soil module TERRA_ML (see p. 111 of Doms et al., 2011) as follows: where T r the transpiration from vegetation, E b and E i are, respectively, the evaporation rates from the (water-permeable) bare-soil fraction f b and that from the intercepted water on the bare soil and vegetation and E snow the snow evaporation.The plant transpiration follows that of the biosphereatmosphere transfer scheme (BATS) of Dickinson (1984).It accounts for both the resistance of water vapour transport from the foliage to the canopy air and the resistance for water vapour transport from the canopy air to the air above the canopy.More details can be found in Ch. 10.2.2 of Doms et al. (2011).
In addition to the standard surface module TERRA_ML, the evaporation from the water storage E imp on the impervious surface fraction Herein, E b is calculated as follows: where E p is the potential evaporation (when positive) or condensation (when negative), H is the Heaviside function which yields 1 in case of evaporation (E p ≥ 0) and 0 in case of condensation (E p < 0), and F m the maximum moisture flux that the soil can sustain (Dickinson, 1984).The potential evaporation E p is obtained from where r H SA is the surface-layer transfer resistance for scalars (heat and moisture), q v and u a are, respectively, the specific humidity and absolute wind speed at the half level of the first model layer and q sat (T s ) is the saturated specific humidity at the surface temperature T s .The calculation of E imp is covered in Appendix A3.The evaporation from the water storage E imp on the impervious surface fraction is calculated as follows: The evaporative area fraction δ imp of the impervious area fraction f imp is calculated according to Wouters et al. (2015): where δ m is the maximum evaporative surface fraction and w m the maximum water storage capacity on the impervious surface.Equation (A10) was obtained by assuming a probability density function of water reservoirs for which the density constantly decreases with increasing reservoir depth.
The amount of water w on the impervious surface changes as follows: where P r is the rain rate, C = −E p H (−E p ) is the condensation, R imp,runoff is the runoff to sewerage and rivers and R imp,infil is the infiltration rate of water into the natural soil originating from the impervious surface (for example, by means of infiltration tubes).The runoff and/or soil infiltration of water hitting the impervious surface occurs as soon as the maximum water storage w m is exceeded: where c imp,runoff is defined as the runoff index of the impervious surface, for which the value depends on the presence of infiltration systems.A value equal to one means that all rainwater exceeding the maximum water-storage threshold from streets and roofs w m is directed to the sewerage and rivers, whereas a value equal to zero means that it infiltrates into the ground.
Default values for w m = 1.31 kg m −2 and δ m = 0.12 are taken from Wouters et al. (2015).By default, all water overshoot is considered as runoff, hence c imp,runoff is set to one by default, meaning that all rainwater overshoot from streets and roofs is considered to be directed to sewerage and rivers.

A4 Anthropogenic heat emission
TERRA_URB accounts for AHE from human activity, which includes the energy dissipation from combustion and electricity consumption.It originates from heating and cooling (such as air conditioning) of buildings and traffic, but also from domestic, industrial and agricultural activity (Sailor, 2011).The AHE is calculated following the methodology of Flanner (2009); see also http://www.cgd.ucar.edu/tss/ahf/.It takes into account latitudinally dependent seasonal and diurnal distribution functions (Flanner, 2009, their Figs. 1 and 2) that are superimposed on the annual-mean anthropogenic heat flux (AHF).The AHE is added to the surface turbulent heat release to the atmosphere, hence it acts as an additional heat source to the first atmospheric model layer.

A5 Tile approach
TERRA_URB makes a distinction between the urban canopy and natural land cover.This is done for each grid cell with a tile approach.Herein TERRA_URB is called twice, once for the urban canopy and once for the natural land cover specifying a different set of bulk parameters.On the one hand, the surface input parameters for the urban canopy are obtained from urban canopy parameters with SURY.On the other hand, surface input parameters for the natural land cover are provided by standard input parameters of the COSMO(-CLM) model.In this way, the ground heat and moisture transport and land-atmosphere exchanges in terms turbulent transport of momentum, heat and moisture are determined separately for each tile.The coupling to the atmospheric model is achieved by weighting each of the land-atmosphere fluxes according to the fractions of the urban canopy and natural land cover.The radiation exchanges are determined by grid cell averaged value of albedo and emissivity, which is weighted according to the respective tile fraction.

A6 Surface input parameters
In addition to the land-surface parameters for the standard COSMO(-CLM) model without urban land-surface parametrization described in Sect.2.2, the implementation of TERRA_URB requires the two additional fields.They include the total ISA and annual-mean AHF.
The ISA field is obtained from the European soil sealing data set representative for the year 2006 of the European Environmental Agency at 100 m resolution (Maucha et al., 2010).Optionally for regions outside Europe, the NOAA data set (Elvidge et al., 2007), globally available at 1 km resolution, can be selected as well.The latter offers flexibility for the extension with upcoming improved global products.In the default configuration of TERRA_URB v2.0, f imp is equal to one for the urban canopy tile (all impervious) and zero for the natural land cover tile (none impervious).As a result, the fraction of urban canopy at each grid cell is equal to the value of ISA and the fraction of the natural land cover is its complement.All vegetation is resolved inside the natural land cover tile, hence the vegetation abundance from the standard input parameters is increased for the natural land cover tile according to its fraction and set to zero for the ur-ban canopy tile.All other parameter values specified for the natural land in TERRA_URB are in accordance to the original version of TERRA_ML and the COSMO(-CLM) model.As a result, the processes (regarding the ground heat and water transport and land-atmosphere interactions) resolved for the natural land cover have not been altered.The only exception is that the vertical profile formulation of the bulk heat capacity (Eqs.5-6) and heat capacity (Eqs.7-8) from SURY is also employed for the natural tile.This was done for providing consistency between the urban canopy tile and the natural land cover tile.Therefore, the height of the natural roughness elements replaces the building height parameter.Default values for the roughness-element height and SAI for the natural land are equal to 0.01 m and 2, respectively.Such a formulation also provides consistency with the TKE-based surface-layer turbulent transfer scheme in which SAI is also equal to 2.
The AHF field is obtained from the global data set of Flanner (2009), which is generated from country-specific data of energy consumption from non-renewable sources.This was apportioned according to population density (conserv-ing the national total)and converted to annual-mean gridded energy flux at a resolution of 2.5 × 2.5 min.By default, a data set is used for which these annual-mean values on the model grid are redistributed according to the ISA field at a scale of 50 km.Therefore, it is assumed that areas with large ISA fraction (including industrial areas with low population densities) have higher anthropogenic heat emissions.For Brussels, the Belgian capital in the centre of the domain in Fig. 1, the AHE reaches an annual-mean value of 49.16 W m −2 in the city centre.These values are of the same magnitude to those obtained in Van Weverberg et al. (2008), i.e. 43.6 W m −2 .While not considering this as a formal validation, the similarity of magnitude of the results obtained by two very different methods inspires confidence.According to the latitudinally dependent seasonal and diurnal distribution functions adopted from Flanner (2009), the values for summer (July) vary between 18.40 W m −2 (night-time) and 37.7 W m −2 (morning and afternoon).For winter (January), they vary between 41.1 and 84.1 W m −2 .
The Supplement related to this article is available online at doi:10.5194/gmd-9-3027-2016-supplement.ceived funding from the Belgian Science Policy Office through its Science for a Sustainable Development Programme under contracts SD/CS/041/MACCBET and BR/143/A2/CORDEX.be.It has also received funding from the Flemish regional government through a contract as an FWO (Fund for Scientific Research) post-doctoral position.The computational resources and services used in this work were also provided by the Hercules Foundation and the Flemish Government (Department of Economics, Sciences and Innovation -EWI).Additional assistance of resources was provided at the NCI National Facility systems at the Australian National University through the National Computational Merit Allocation Scheme supported by the Australian Government.We would like to thank Gianluca Mussetti (Swiss Federal Laboratories for Materials Science and Technology -EMPA) and Johan Zürger (Austrian Institute of Technology -AIT) for testing TERRA_URB for regional-climate modelling over Zürich and Vienna.We also acknowledge Jürgen Helmert and Burkhardt Rockel for providing support in extending the EXTPAR input parameter tool, and the COSMO consortium and CLM community for their support.Finally, we wish to thank Chunlei Meng, Yi-Ying Chen and the other reviewers for their insightful suggestions to the manuscript.
Edited by: M.-H.Lo Reviewed by: C. Meng, Y.-Y.Chen, and three anonymous referees

Figure 1 .
Figure 1.Domain composite of the impervious surface area, vegetation cover and orography (indicated with the shading) for the reference (REF) model setup of the COSMO-CLM model coupled to TERRA_URB v2.0 at 2.8 km resolution over Belgium.The arrows directed upwards indicate the locations of the in situ observations located in the urban area of Antwerp (left; Royal Lyceum of Antwerp) and in the rural area of Vremde (right; Organic Farm van Leemputten).The arrows directed downwards indicate the locations of the tower observations located at a flat industrial terrain in Zwijndrecht (left arrow) of the Flemish Environmental Agency and at an rural area (right arrow) of the Belgian Nuclear Research Center (SCK•CEN) in Mol.

Figure 2 .
Figure 2. Modelled horizontal profiles of screen-level (top panels) and surface (bottom panels) temperatures at noon (left panels) and midnight (right panels), averaged for 21 July to 20 August 2012.

Figure 3 .
Figure 3. 24 h running averages and mean diurnal cycle for modelled (thick lines) and observed (OBS; blue stars) temperatures during mid-summer (21 July to 20 August 2012) at the Royal Lyceum for Antwerp (urban), the Organic Farm Van Leemputten (rural) and their difference (the urban heat island effect of Antwerp).The simulation with the COSMO-CLM model coupled to the standard land-surface module TERRA_ML without urban parametrization is indicated with STD (black), whereas the reference simulation with the COSMO-CLM model plus TERRA_URB v2.0 with urban parametrization is indicated with REF (red).The dotted lines indicate the range between the 16th and 84th percentile of the observed temperatures, whereas the grey and light red areas indicate the ranges for the simulations STD and REF, respectively.An overview of the canopy parameter sensitivity simulations (AL, AH, BL, BH, etc.) can be found in Table2.

Figure 4 .
Figure 4. Evaluation of modelled land-surface temperatures against MODIS satellite observations.A distinction is made between four urban classes: no urban (NU) with impervious surface area (ISA) ≤ 0.05; light urban (LU) with 0.05 < ISA ≤ 0.25; medium urban (MU) with 0.25 < ISA ≤ 0.5; dense urban (DU) with ISA > 0.5.The amount of valid pixels are shown above each day and night plot.The violin plots represent the distribution of LST for MODIS in blue and the reference simulation (REF) with COSMO-CLM coupled to TERRA_URB v2.0 in orange.Full and dashed lines represent the median for the low (AL-GL) and high (AH-GH) scenarios from Table2, respectively.The upper panels show the sample distributions of the absolute temperatures, whereas the lower panels show the sample distribution of the difference between median of each urban class (LU, MU and DU) with the NU class.

Figure 5 .
Figure 5. Observed (stars) and modelled (lines) nocturnal (0 h) vertical profiles for VMM industrial site in Zwijndrecht and the SCK•CEN rural site in Mol, averaged for the summer period of 21 July to 20 August 2012.STD (black) indicates the simulation with standard model COSMO-CLM without urban parametrization.REF uses the same setup as STD, but showing results of TERRA_URB v2.0 coupled to COSMO-CLM.An overview of the remaining simulations (AL, AH, BL, BH, etc.) can be found in Table 2.The dotted lines indicate the range between the 16 and 84th percentile of the observed profiles, whereas the grey and light red areas indicate the ranges for the simulations STD and REF, respectively.

Figure 6 .
Figure 6.Application outline of the SURY.The main application flow of SURY is indicated with full arrows, whereas optional application flows are indicated with dashed arrows.

Table 1 .
Loridan and Grimmond (2012)rban canopy parameters.They are taken as input for the SURY.The default urban canopy parameters correspond to the recommended values for the medium urban category inLoridan and Grimmond (2012); see their Table4(stage 5b).The lower panel shows the bulk parameters, which is the output of SURY.The parameter u * refers to the friction velocity.

Table 2 .
Overview of parameter sensitivity experiments.Seven couples of experiments (AL, AH, BL, BH, CL, CH, DL, DH, EL, EH, FL, FH, FL, GH) are performed for which the default urban canopy parameters are modified to the values in the low (L) and high (H) columns.
as FL09 in the table) are used.

Table 3 .
Daily averages for modelled and observed temperatures (unit: K) during mid-summer (21 July to 20 August 2012) at the Royal Lyceum for Antwerp (urban), the Organic Farm Van Leemputten (rural) and their difference (the urban heat island effect of Antwerp).Each first row shows results of the observations.Each second row shows results for COSMO-CLM model without urban parametrization (L column), for COSMO-CLM plus TERRA_URB v2.0 (H column) and their absolute difference (| H − L | column).The remainder rows show each of the urban-parameter sensitivity simulations in Table2.Therefore, the low scenarios are shown in the L column, the high scenarios in the H column and the parameter sensitivity range (absolute difference between L and H) in the| H − L | column.

Table 6 .
Modelled and observed temperatures (in K) vertically averaged for the nocturnal (0 h) tower profiles at the VMM industrial site in Zwijndrecht and at the SCK•CEN rural site in Mol.− L | column).The remainder rows show each of the urbanparameter sensitivity simulations in Table2.Therefore, the low scenarios are shown in the L column, the high scenarios in the H column and the parameter sensitivity range (absolute difference between L and H) in the | H − L | column.