The role of urban boundary layer investigated by high resolution models and ground based observations in Rome area : a step for understanding parameterizations potentialities

abstr The urban forcing on thermo-dynamical conditions can largely influences local evolution of the atmospheric boundary layer. Urban heat storage can produce noteworthy mesoscale perturbations of the lower atmosphere. The new generations of high-resolution numerical weather prediction models (NWP) is nowadays largely applied also to urban areas. It is therefore critical to reproduce correctly the urban forcing which turns in variations of wind, temperature and water vapor content of the planetary boundary layer (PBL). WRF-ARW, a new model generation, has been used to reproduce the circulation in the urban area of Rome. A sensitivity study is performed using different PBL and surface schemes. The significant role of the surface forcing in the PBL evolution has been verified by comparing model results with observations coming from many instruments (LiDAR, SODAR, sonic anemometer and surface stations). The crucial role of a correct urban representation has been demonstrated by testing the impact of different urban canopy models (UCM) on the forecast. Only one of three meteorological events studied will be presented, chosen as statistically relevant for the area of interest. The WRF-ARW model shows a tendency to overestimate vertical transmission of horizontal momentum from upper levels to low atmosphere, that is partially corrected by local PBL scheme coupled with an advanced UCM. Depending on background meteorological scenario, WRF-ARW shows an opposite behavior in correctly representing canopy layer and upper levels when local and non local PBL are compared. Moreover a tendency of the model in largely underestimating vertical motions has been verified.


Introduction
Numerical weather prediction (NWP) models can, nowadays, work at very high resolutions (∼ km), but there are many subgrid processes that develop at finer scales and can not be explicitly represented.They have to be included into the models for correctly reproducing the atmospheric state.Turbulent mixing is one of the subgrid phenomena that has a large impact on the state of the atmosphere; it occurs in low layers of the atmosphere and it characterizes the planetary boundary layer (PBL), taking charge of the vertical transport of mass, heat and momentum.The relatively high frequency of occurrence of turbulence near the ground differentiates the PBL from the rest of the atmosphere (Stull, 1988).NWP models have to reproduce turbulence at various scales and so they need appropriate representations of the PBL.Many different PBL schemes are available; they differ from each other by the vertical mixing formulation and the closure order.Some parameterizations have computational advantages (like the ones based on the so called local-K approach), but they can fail to reproduce the mass and momentum transport accomplished by large eddies (Stull, 1993).In order to overcome these deficiencies parameterizations accounting for nonlocal contributions by countergradient terms (nonlocal-K approach) or more sophisticated representations like higher-order closure approaches based on prognostic prediction of turbulent kinetic energy (TKE) have been developed.Previous studies have demonstrated the important role of PBL schemes for improving models prediction skills; E. Pichelli et al.: Planetary boundary layer of the urban area of Rome they can have an impact of precipitation forecast (Hong and Pan, 1996) and influence near-surface wind representation (Grossman-Clarke et al., 2008).The accuracy of the PBL schemes affects forecasts of both local and large-scale meteorological phenomena (Hacker and Snyder, 2005).It has been demonstrated, moreover, that the representation of characteristics of the PBL is more or less sensitive to different parameterizations, depending on whether one considers the mean or turbulent structure of the layer (Holt and Raman, 1988).
Recently, the new generations of high-resolution NWP models have also been applied to urban areas for both weather forecast and research purposes (Grossman-Clarke et al., 2008;Salamanca et al., 2011Salamanca et al., , 2012;;Flagg and Taylor, 2011;Kim et al., 2013;Wouters et al., 2013).Generally, these areas are covered by dry materials (asphalt or concrete) and are warmer and drier than adjacent rural areas (Oke, 1982); the efficiency of the urban areas in storing heat can produce remarkable mesoscale perturbations of the lower layers of the atmosphere (by variations of wind, temperature and water vapor content).The inclusion of urban effects in NWP models can, thus, significantly improve the weather forecast (Kim et al., 2013).Featuring the urban boundary layer adds a further complexity that is usually resolved by coupling PBL schemes with urban canopy parameterizations inside models.Collier (2006) clearly stated the need for a better understanding of the PBL of the urban areas because of their impact on the weather.
The purpose of this study is twofold: 1. to investigate the local circulation in the urban area of Rome using both observations from different instruments and the new-generation WRF (weather research and forecasting) model (Skamarock et al., 2008) simulations; 2. to investigate the model's ability to reproduce the circulation in the urban area of Rome.
To the authors' knowledge ,this study is the first one to analyze the horizontal and vertical structure of the urban area of Rome using a set of measurements and the high resolution model WRF.The position of this city very close to the Tyrrhenian coast (west side of italian peninsula) make the local atmospheric circulation complex because of the interactions between urban heat island (UHI) and sea breeze; observational works (Colacino and Dell'osso, 1978;Mastrantonio et al., 1994;Casadio et al., 1996) have established the importance of the sea breeze in the Rome area.A few previous numerical studies (Ferretti et al., 2003;Caballero and Lavagnini, 2002) show that simulation of local circulation in Rome area is challenging for the model.The approach of this study is based on an the extensive literature already published that investigates the response of the PBL schemes available for models (MM5 at first and then WRF), with respect to different meteorological events on urban or rural areas (Dandou et al., 2005;Grossman-Clarke et al., 2008;Thomsen and Smith, 2008;Trusilova et al., 2008, etc.).Previous studies allowed for both highlighting biases of the model and for understanding the errors mechanisms' generation (Dandou et al., 2005).Based on our experience (operational run) WRF responds correctly in cases of strong forcing, but in summertime it may erroneously represent local circulation when large-scale influence is weak.The WRF model offers numerous options for PBL schemes; recent studies (Shin and Hong, 2011) have compared some of them, concluding that nonlocal schemes are more favorable under unstable conditions, whereas TKE closure schemes perform better under stable conditions, even with a large bias for most variables.They show that main differences that arise are due to the local or nonlocal nature of the parameterizations; the one exception is for near-surface variables, which, on the other hand, are strongly influenced by surface schemes, which in turn significantly influence the structure of the PBL (Pan and Mahrt, 1987;Stull, 1988).
In the present work WRF has been tested for the Rome area using local and nonlocal schemes.Among the ones available for WRF, we have chosen the Yonsei University (YSU, Hong et al., 2006) scheme as the nonlocal scheme, it being the new generation of the Medium Range Forecast (MRF, Hong and Pan, 1996) of MM5 model, largely used both in operational simulations for the area of interest and for specific studies (Ferretti et al., 2003) with good results.On the other hand, Shin and Hong (2011) shows that local schemes tend to converge to similar results for turbulent structure of the PBL and that surface schemes are responsible for differences near the surface more than PBL scheme.The Mellor-Yamada-Janjic (MYJ, Mellor and Yamada, 1982) scheme has been chosen as TKE closure scheme; this is one of the most-used schemes and it can also be coupled with the multilayer canopy model available with WRF (Martilli et al., 2002).It has indeed been shown (Flagg and Taylor, 2011;Lee et al., 2011;Kim et al., 2013;Martilli, 2002) that an appropriate explicit parametrization of urban physical processes produces more accurate results in other urban areas.Based on the previous considerations, the sensitivity to the surface layer of the two PBLs is also investigated in order to highlight their role in influencing the PBL inside and outside the Rome urban area.As a further step, simulations with different UCMs have been performed.
Observations from sonic anemometer, lidar, sodar, and ground-based stations have been used for the comparison presented.A brief description of instruments, experimental techniques, and the model configuration is given in the first part of the paper.An analysis of the case study and the comparison between the model and the observations is presented in the second part.These allow for an investigation of the circulation in the Rome area in a meteorological scenario characterized initially by local forcings and weak large-scale influence and later by weak convection conditions due to a moderate wind regime driven by large-scale forcing.A brief comparison between the model and the observations in rural sites is also presented to investigate possible different impacts of PBL, land-surface models (LSMs) and UCMs outside the urban area.

The Sonic anemometry
Sonic anemometry is mainly used in atmospheric turbulence studies.This instrument allows one to measure threedimensional wind velocity and sonic temperature.From the latter it is possible to retrieve the virtual air temperature.The operating principle is that the time lag of the sonic waves propagating in moving air depends on air speed and direction; it measures the "transit time", i.e., the time it takes for ultrasonic signal to travel from one transducer to another.An ultrasonic anemometer (model 81 000 V of the Young Company, USA) was installed in 2007 on the roof of the building of physics department at the Sapienza, University of Rome (41.9 • N, 12.5 • E, 75 m a.s.l.).This place is located in the city center which is a very populated area, strongly influenced by anthropogenic activity (Meloni et al., 2000).Wind speed is measured with an accuracy of ±1 % in the range of velocities 0-30 m s −1 and of ±3 % in the range between 30 and 40 m s −1 .Wind direction accuracy ranges between ±2 • (0-30 m s −1 ) and ±5 • (30-40 m s −1 ).The air temperature (T a ) and relative humidity (RH) probes combine thermistors (accuracy: ±0.15 • C; range of measure: −30/50 • C; response time: 20 s) and a capacitive hygrometer (accuracy: 3 % in the RH range of 10-100 %).Both of the meteorological sensors and the sonic anemometer are connected to a data logger.The rate at which the T a and RH are sampled via anemometer is 32 Hz.Horizontal and vertical wind components together with sonic temperature and the RH and T a values are sampled every 30 min.Spikes were identified and removed when the readings were above/below ±3.5 std (Hojstrup, 1993;Vickers and Mahrt, 1997) taking into account a temporal window of 100 s.Reprocessing software performs a quality control of the data set and provides every 30 min the means of three wind components and their standard deviations, the horizontal mean wind, the mean direction, the mean sonic temperature and its standard deviation.The same for T a and RH.In addition, some turbulence parameters, such as friction velocity and turbulent heat transfer, are retrieved.

The lidar system
The system has been designed to observe atmospheric aerosol vertical profiles in the PBL and the free troposphere.The radiation source is a Q-switched single-stage Nd : YAG (neodymium : yttrium-aluminum-garnet) laser emitting linearly polarized pulses at 1064 and 532 nm, with a repetition rate of 10 Hz.The 532 nm radiation is produced by a second harmonic generator crystal.The backscattered radiation is collected by a 100 mm diameter reflector telescope in a Cassegrain configuration (Thomas, 1995).The laser beam is directed towards the zenith coaxially to the Cassegrain receiver, hence the telescope secondary mirror masks the strong atmospheric echoes from the lower 300 m; in this way, saturation of the detectors is prevented and the receiver sensitivity and field of view (FOV) is regulated to observe the atmosphere up to the tropopause.Another small-aperture, large-FOV refractor telescope receiver is placed beside the Cassegrain, but sufficiently close to the laser beam to observe the strong echo from the lowest atmosphere.The collimated signals are filtered by narrow-band interference filters to reduce the sky light and to allow measurements even in full daylight during the summer.Simultaneous analog detection and single-photon counting is performed on the signals received by the larger telescope, while analog detection is only performed on the signal from the smaller receiver.The low-range analog signal and the high-range analog and photon counting signals are matched in the overlapping altitude ranges to produce a continuous lidar signal from about 50 m to the upper troposphere with a vertical resolution of 7.5 m.The acquisition system is programmed to perform an integration of the backscattered signals over 300 laser shots (corresponding to 30 s).This is the highest time resolution achievable from the saved raw data but all the analyses were performed on profiles averaged over 5 min (corresponding to 3000 laser shots).The retrieval of the backscatter ratio R, defined as the ratio between the total (aerosol and molecules) backscattered signal, and the portion due to the atmospheric molecules only, is obtained by a standard algorithm.The procedure for estimating the height of PBL (PBLH) is based in the computation of three parameters vs. time as follows.
1. Three derived quantities from the profiles of R(t, z), namely TV (time variance), FD (first-order derivative in z) and SD (second-order derivative in z) are calculated using a three-point Lagrangian interpolation along the vertical coordinate.Both FD and SD are smoothed by a running average over 75 m in the vertical (this is the error associated with PBLH.The procedure to achieve the final SD profile requires two smoothing operations).
2. From the closest-to-the-ground relative maximum in the TV vertical profile a negative maximum in the FD followed by a positive maximum in the SD are searched for in a range of 400 m of height.If the search is successful, the three heights for TV, FD and SD are saved (TVH, FDH and SDH respectively)

The sodar system
The Doppler sodar is operated in a three-axis monostatic mode with a pulse repetition period of 6 s, allowing a maximum probing range of 1000 m.Two antennae are tilted 20  (Argentini et al., 1992).The sodar detects the vertical profiles of the three components of the wind every 60 s and vertical profiles of the turbulence intensity every 6 s.A data processing similar to the one used for the lidar is applied to the turbulence intensity profiles in order to retrieve the height of the mixing layer height.

The Model configuration
The nonhydrostatic WRF-ARW model (Advanced Weather Research and Forecasting) V3.4.1 is used for this study; this is a primitive equations model with a terrain following vertical coordinate and multiple nesting capabilities (Skamarock et al., 2008).Four two-way nested domains are used (Fig. 1) to enhance the resolution over the urban area of Rome and its surroundings.The mother domain is centered at 41.116 -35 unequally spaced vertical levels, from the surface up to 100 hPa, with a higher resolution in the planetary boundary layer than in the free atmosphere and the lowest level at a 47 m of height; -long-wave RRTM (Rapid Radiative Transfer Model) and short-wave Dudhia schemes for radiative transfer processes.Both these parameterizations come from the MM5 model and are respectively based on Mlawer and Dudhia-Lacis-Hansen schemes; -Kain-Fritsch cumulus convection parameterization is applied to domains 1 and 2; whereas no cumulus scheme is used for domains 3 and 4; -Morrison two-moment bulk scheme for microphysics; -US Geological Survey (USGS) land use.
Numerical experiments have been performed using different PBL parameterizations and different combinations of PBL schemes and surface models with the aim of both assessing the correct configuration for the urban area of Rome and investigating its local circulation.It has to be pointed out that the urban area of Rome has not experienced large transformations in the last century, whereas the suburban area has considerably changed.The USGS land use is not the most upgraded data set, but it is still a good choice for the representation of the urban area where the study is focused.
The following parameterizations are used for boundary layer: -the Yonsei University (YSU) scheme (Hong et al., 2006), based on the nonlocal K-theory mixing in the convective PBL (Troen and Mahrt, 1986); -the Mellor-Yamada-Janjic (MYJ) PBL (Mellor and Yamada, 1982;Janjic, 2002); this is a local 2.5 turbulence closure model, with an upper limit imposed on the length scale that depends on the turbulent kinetic energy (TKE) as well as on the buoyancy and shear of the driving flow.
To account for the surface physics, a surface scheme together with a land-surface parameterization is needed.The first computes friction velocities and exchange coefficients that enable the calculation of surface heat and moisture fluxes by the land-surface models.These fluxes provide a lower boundary condition for the vertical transport in the PBL; the land-surface model update the land's state variables.Two different schemes are used for both surfaces (Skamarock et al., 2008): -the MM5 surface model based on similarity theory (MO-MM5); -the Eta surface layer (MOY-MYJ); and land surface: -the MM5 five-layer thermal diffusion scheme (TD-MM5); -the Noah land-surface model (NoahLSM).
With the aim of representing the city scale effects at the mesoscale, the NoahLSM is also coupled to a UCM; a singlelayer scheme (UCM1) is coupled with YSU PBL (Kusaka et al., 2001), whereas both a single-layer (UCM1) and a multilayer model (UCM2) are used with MYJ scheme (Martilli et al., 2002).As the MYJ simulation with UCM1 does not produce remarkable differences with respect to the one without any urban canopy scheme coupling, only the simulation with the UCM2 will be presented in this paper.Only one category is used to characterize urban land use; default urban morphology parameters (street width, building heights, etc.) associated with the urban category have been considered representative, on average, of the Rome environment.These are in line with settings of Kim et al. (2013) for the city of Paris.The largest difference is found for roof width and anthropogenic heat flux (AH); the first one, set by default at 9.4 m (3.75 m for Paris), is close to a realistic average between the roof of the building where instruments site is located and surroundings ones.The maximum AH is left at 50 W m −2 (Paris ranges between 35 and 70 W m −2 in May).As a future step the authors will further "adjust" the urban parameters, but the configuration here used can be considered representative of Rome area even if not totally realistic.
In order to highlight the sensitivity of the nonlocal PBL scheme with respect to the land-surface several simulations using the YSU PBL parameterization are performed: (1) using the TD-MM5 for the land-surface (YSUtd); (2) using Noah land-surface model (YSUNoahNOURB); (3) using the same configuration as option 2, but adding the urban canopy model (YSUNoahUCM1).The same set of simulations are performed using the local 2.5 turbulence closure MYJ model: (1) the MYJtd using the TD-MM5 for the land-surface; (2) the MYJNoahNOURB using Noah land-surface model; (3) using the same configuration as option 2, but adding urban canopy model (MYJNoahUCM2).In Table 1 a summary of simulations performed with different configurations is shown.The acronyms in the first column will identify each simulation hereafter.The ECMWF (European Centre The following meteorological parameters are used for this analysis: measurements at 25 m of the horizontal wind velocity (m s −1 ) and direction (deg) as well as vertical velocity (m s −1 ), friction velocity (m s −1 ) by the sonic anemometer and temperature ( • C) and relative humidity (%) by combined probes; the PBL height (m) time series retrieved by the lidar measurements; the time series of the horizontal and vertical wind profiles measured by the sodar.It has to be pointed out that the anemometer detects the atmosphere within the canopy layer, whereas the sodar and lidar both scan the atmosphere above this layer, but within the PBL.Therefore, the anemometer measurements are not directly comparable with the sodar and lidar ones, but coupling the two allows for investigating a large part of the PBL.

Observations
The 6-7 February 2008 case is characterized by the transition from a weak wind regime to one of weak convection associated to moderate wind.During 6 February, a low-pressure Table 2. Anemometer mean standard deviations for 6-7 February 2008; wsp indicates horizontal wind speed, wdr is the horizontal wind direction, w is the vertical wind velocity, u * is the friction velocity, T is the temperature, RH is the relative humidity.6-7 February 2008 wsp (m s −1 ) wdr (deg) w (m s −1 ) u * (m s −1 ) T ( • C) RH (%) STD 0.9 29 0.5 0.06 0.2 0.9 system in the southeast of Italy is moving southeastward and an anticyclone is entering from Gibraltar moving northeastward, thus producing a weak wind regime over central Italy (Fig. 2a).The satellite images (not shown) reveal clear sky all day over the urban area of Rome except for some cloudiness at night.The following day (7 February), the anticyclone reaches central Europe and the cyclone extends over the south Mediterranean; this produces an increase in the surface pressure gradient and of the easterly wind at the middlelow atmosphere over Italy (Fig. 2b).This is associated with an outbreak of cold and dry air, which triggers weak but nonprecipitating convection, mixing the lower atmosphere.Clear sky also characterizes the urban area of Rome on 7 February (not shown).
The horizontal wind component is weak, as shown by the sonic anemometer during the first day, except for a maximum (Fig. 3a, red line) recorded at approximately 15:00 UTC.An increase of the wind strength is observed before midday of the second day, lasting for the whole period.The large-scale structure would suggest that the two maxima are produced by different mechanisms: the first one by a thermally driven circulation, the second one by the large-scale forcing.
The wind direction (Fig. 3b, red line) shows a westerly/northwesterly wind during the second part of 6 February, probably produced by the interaction of large-scale flow with the onset of the sea breeze as the timing and the persistence (from 12:00 to 18:00 UTC) of this wind regime would suggest.Sea breeze interaction with the circulation in the urban area of Rome is a well-established phenomenon, even in wintertime (Mastrantonio et al., 1994;Ferretti et al., 2003).During the second day a larger variability is detected; this is associated with a north westerly wind (between early morning and 01:00 UTC), in phase with an updraft registered during 7 February (Fig. 3c, red line) and an increase of the frictional velocity (Fig. 3d, red line).
Concerning vertical velocity, it has to be pointed out that anemometer measurements are affected by large errors; standard deviation (Fig. 3c, gray bars) shows values comparable with the measurements itself (∼0.5 m s −1 ).The mean standard deviations time series for 6-7 February for each anemometer variable are presented in Table 2. Two diurnal updrafts are detected by the anemometer due to both free and forced convection (Fig. 3c, red line): at approximately 15:00 UTC on 6 February and at 11:00 UTC of the day after.The first updraft (Fig. 3c, red line at 15:00 UTC on 6 February) is mainly due to a plume which is sustained by the thermal contribution of the UHI and also by mechanical contribution, as suggested by the friction velocity increase well in phase with this updraft (Fig. 3d, red line).During the second day, larger velocities are detected for both the horizontal and vertical wind component (Fig. 3a and c, red line, after 11:00 UTC on 7 February) because of an increase of the instability produced by the cold and dry air outbreak; moreover, a second updraft is registered after 18:00 UTC of comparable intensity with the diurnal one.This allows for the inference that, during the second day, the increase of the horizontal wind partially inhibits free convection, but contributes to the development of local mechanical turbulence.The PBL height retrieved by lidar supports the previous hypothesis.It clearly shows a well-developed PBL during the first day (Fig. 7a, red line), supporting the hypothesis of meteorological conditions strongly driven by the local forcing, that is by a UHI.On the other hand, the large variability and the shallower PBL recorded during the second day (Fig. 7a) confirms an inhibition of free convection due to the large-scale circulation preponderance.The time series of the horizontal and vertical wind vertical profile detected by sodar also support the previous hypothesis.The vertical profile of the horizontal wind speed (Fig. 8a) shows a weak signal during 6 February, between 50 and 300 m of height.After midnight a wind increase is recorded, with a maximum greater than 8 m s −1 above 100 m and strong winds reaching 14 m s −1 during the second part of 7 February.On the other hand, the vertical wind profile (Fig. 10a) shows positive values for the whole period, except for a downdraft after 12:00 UTC on 6 February, whereas an intermittent signal is recorded by midday of 7 February.Values ranging between 0 and 80 cm s −1 are detected by sodar (Fig. 10a), high variability is also found during 6 February: two main updrafts develop between 10:00 and 24:00 UTC, with maxima at 11:00 and 18:00 UTC; the updrafts are mostly at upper levels, and are associated with two weak and short downdrafts after midday.During 7 February mainly upward motion is detected with two maxima, one between 01:00 and 05:00 UTC, and one between 10:00 and 16:00 UTC.

WRF results
As a first step, the model results are compared with the local observations detected in the urban area for 6-7 February 2008.The meteorological parameters detected by the sonic anemometer and connected probes are compared with the one produced by WRF interpolated at the same height (25 m).The PBL height time series retrieved by lidar measurements and the time series of the horizontal and vertical wind profiles detected by sodar are compared with the one produced by WRF at the same location.All the model results are analyzed at the highest resolution (0.78 km).

Inside the urban area
The comparison between WRF and the anemometer (Figs.YSU during the first hours of 6 February except for a well reproduced maximum after noon.During 7 February, WRF is able to reproduce the wind regime change, but a large error is found during the phase in which the wind speed increases; this turns into an overall positive bias of the model during the moderate winds regime.During this phase the YSUtd (blue line) and the YSUNoahNOURB (yellow line) tend to produce similar results, whereas the YSU coupled with the urban canopy scheme (YSUNoahUCM1, Fig. 3a, green line) produces the largest errors (≈ 8 m s −1 ).The WRF overestimation of horizontal wind can be produced by an excess of forcing of the upper layers' dynamic to the urban canopy layer, as the comparison with sodar vertical profiles would suggest (Figs. 8 and 9).As will be pointed out later, the YSU wind speed maxima are developed at lower levels than the sodar ones, strongly supporting the previous hypothesis.
The wind direction (Fig. 3b) is partially well-reproduced by the YSU simulations: the model correctly reproduces and maintains the sharp change of direction registered after midday on 6 February and after 18:00 UTC on 7 February.Errors are found during the remainder of the two days, especially in the early morning and late afternoon of 6 February.Biases between YSU and the anemometer range between ∼ 36 and ∼ 177 • , with the largest error associated with the YSUNoahUCM1.During nighttime (from 22:00 UTC on 6 February to the early morning of 7 February), the anemometer measurements (Fig. 3b, red line) show mainly a northwesterly wind, whereas the model produces wind initially coming from northeast (reproducing the typically night flow of the Tiber valley); after 06:00 UTC on 7 February, a sharp change in the wind direction in agreement with measurements is produced (Fig. 3b, blue, yellow and green lines).Also for the wind direction, no relevant differences are found using UCM1 (Fig. 3b, green line).
Concerning the vertical velocities, the comparison between the YSU simulations and the anemometer (Fig. 3c) shows a large underestimation of the mean vertical winds at low levels during most of the simulation, regardless of the surface scheme (Fig. 3c, blue, yellow and green lines).With the aim of understanding the model ability in reproducing local forcing, the WRF vertical velocity is multiplied by a factor of 10 in order for the researcher to be able to compare the observed and modeled signals.Besides the large underestimation, all YSU simulations produce a downdraft instead of an updraft during daytime of 6 February, because of an anticipation of the midday plume with respect to observations.On the other hand, the midday maximum on 7 February is well in phase with the observed one, whereas the late evening one is delayed.This would suggest a model difficulty in reproducing the different contributions (mechanical and thermal) to the updraft.The friction velocity analysis will help clarify this point.The comparison between the YSU friction velocity and the anemometer (Fig. 3d) shows a fairly good agreement during the first day of weak advection regime, although a time anticipation of the 6 February maximum is produced.On the other hand, although the increase of the frictional velocity reproduced by WRF is well in phase with the observed one, a large bias is found during the second day.
The comparison between the anemometer and the WRF temperature and relative humidity (Fig. 4) shows a tendency to underestimate the maxima of both variables.All WRF simulations correctly reproduce the rate of T a increase (Fig. 4) during daytime for both days, but the maximum on 7 February is largely underestimated regardless of the surface scheme.The attempt to reproduce the temperature change during the night is noteworthy; the anemometer shows a minimum by 01:00 UTC that lasts for approximately 4 hours, decreasing slightly by 06:00 UTC.The model is able to reproduce this signal but it anticipates the minimum by approximately 2 h; on the other hand, it is able to produce the plateau lasting for 4 h in good agreement with the observations, but then produces a negative bias for the 06:00 UTC minimum.The decreasing rate of temperature is fairly well-reproduced for both days, but the timing advance in reproducing daily maximum turns into a temperature underestimation during second part of each day, which is larger for 7 February.The relative humidity (Fig. 4b) is reproduced by WRF with some discrepancy; the model shows poor sensitivity to the land-surface scheme, especially during the first day.The first day maximum is well reproduced by all YSU simulations, with YSUtd producing the best agreement (Fig. 4b, blue line), whereas YSUNoahNOURB and YSUNoahUCM1 slightly underestimate it (Fig. 4b, green and yellow lines) with negligible differences between them.All YSU simulations underestimate the midday minimum for both days and produce them in advance; during the first day this is caused by an overestimation of the RH decreasing rate.During the second day all the simulation reach the right maximum value, but none of them is able to reproduce the plateau between 00:00 and 06:00UTC; they begin the decreasing phase early and thus produce a general underestimation of the observations.
In the following, the comparison between observations and MYJNoahUCM1 is not presented because large similarities with both MYJtd and MYJNoahNOURB are found.The UCM2 is considered instead, as it produces some noteworthy differences regarding the simulation.
The comparison between MYJ simulations and the anemometer  show trends similar to YSU ones, but a better agreement with measurements is found for most variables.During 6 February and early 7 February, all the simulations show a good agreement with the anemometer for the horizontal wind speed (Fig. 5a) producing values within the anemometer error; in addition, the maximum between 12:00 and 18:00 UTC on 6 February is reproduced at the correct time by all simulations, although MYJNoahUCM2 underestimates its extent.On the other hand, an overestimation is found during 7 February for this scheme, regardless of the land-surface scheme, if no urban canopy is activated (Fig. 5a, orange and black lines).A smaller maximum error (6.0 m s −1 ) than YSU one is found.The UCM2 activation (Fig. 5a, pink line) dramatically reduces the horizontal wind intensities, almost canceling the error during 7 February.Urban effects produced by UCM2 thus allows for decoupling low level from those above.
The MYJ is able to reproduce the wind direction time series fairly well during most of the simulation (Fig. 5b), even missing the sharp changes registered by anemometer in the early hours of 6 February, similarly to YSU.During late evening it produces smaller errors than YSU, but a comparable signal is shown at nighttime with wind coming mainly from the east/northeast, whereas the anemometer detected wind from the northwest.
As in the YSU, the vertical velocity is largely underestimated in all MYJ simulations (Fig. 5c) regardless of the land-surface scheme.Also for MYJ the WRF vertical velocity is multiplied by 10 in order to enable one to compare the observed and modeled signals and assess a partial correlation between them.In addition to a time anticipation, an attempt to produce the first-day updraft is shown with increasing values if the UCM2 is activated (Fig. 5c, pink line).Concerning the second day, all MYJ simulations attempt to reproduce the bimodal structure developed during morning and late afternoon.MYJNoahUCM2 further reduces velocities, thus increasing the bias.
The friction velocity for MYJ (Fig. 5d) show trends very similar to YSU ones: a fair agreement is shown during the daily cicle, with a smoother increasing and decreasing rate during 6 February than in YSU, thus allowing for reducing bias with respect to the anemometer.On the other hand, a large disagreement is found during 7 February.In this case some sensitivity to surface model is found; if the urban canopy model is used, better results are obtained for nightime than those obtained with MYJtd and MYJNoahNOURB, but larger errors are shown later.This suggest that a further adjustment of urban parameters should be carried out, but this is left for a future work.
The temperature comparison (Fig. 6a) shows that local parameterization (MYJ) also underestimates observations, regardless of the land-surface scheme, but it reproduces both cooling and warming rates well for both days.If the UCM2 is activated (Fig. 6a, pink line) the errors in the two maxima are smaller than those in the other MYJ simulations, whereas errors comparable to the other two simulations are shown for minima.A good agreement with the anemometer is found also for the relative humidity, but with a larger underestimation of the first day maximum than the YSU (Fig. 6b).The error is recovered after early morning of 6 February, after which MYJ shows an overall better agreement with observations than YSU.The MYJtd (Fig. 6b, orange line) produces the best estimation of RH except in the last part of the simulation, when larger overestimation than that of YSU is found.The previous analysis allows for inferring a WRF tendency to overestimate horizontal wind component at low levels regardless of the PBL parameterization; an excess of interaction with the large-scale structures during strong large-scale forcing conditions is supposed.This would suggest an overestimation of vertical transport of horizontal momentum due to an inefficiency in decoupling the canopy layer from the upper ones.In general, coupling the PBL with Noah LSM does not reduce the model error for the wind.A large reduction of the wind speed error is found for the local PBL scheme if the multilayer urban canopy model is used.The nonlocal scheme, by contrast, shows poor sensitivity to the urban scheme.Both the parameterizations underestimate maximum temperature during daytime at the site height (25 m).The error can be associated with the underestimation of the temperature at lower layers as verified in the comparison with ground-based stations shown in next paragraph and in agreement with findings of Hu et al. (2010) for mean diurnal variation of temperature detected at 2 m above ground level (2 m temperature) in southeast American sites.The local scheme shows greater ability than YSU in reproducing the relative humidity.
The previous results would suggest that the local 2.5-order closure PBL better reproduces the low-levels PBL in urban areas both for a meteorological scenarios characterized by local circulation and by large-scale signal influencing the low levels.In this second case, a multilayer urban canopy scheme allows for reducing the errors for most variables.
To investigate the PBL vertical structure produced by WRF, lidar and sodar measurements are used.Lidar data provides the PBL height, which is a key variable for the parameterizations because it drives the representation of nonlocal mixing.For this case study, PBL height measurements by lidar are available; measurements by sodar are also available for some time intervals (Fig. 7, respectively red and black dashed lines).Lidar measurements (Fig. 7, red line) show a well-developed PBL during 6 February, that reaches a height of 1200 m.During 7 February a lowering of the maximum height is recorded (800 m) associated with a high frequency variability, also suggesting a very turbulent state of the atmosphere during the late afternoon.The large friction velocity values measured by the anemometer support this  hypothesis (Fig. 3d, red line).The PBL height retrieved by sodar (Fig. 7, black dashed line) is available mostly during nighttime and early morning; this is usually more accurate than lidar height measurements below 500 m (errors never exceeding 75 m); for this event a large agreement between the two instruments is found.
The comparison between lidar and WRF shows that the model reproduces most likely the PBL growth during the first day, but underestimates and anticipates (∼ 1.5 h) the maximum if the nonlocal scheme is used (Fig. 7a).Besides the time shift, a bias of about 700 m with respect to observations is found.On the other hand, MYJ (Fig. 7b) reduces the maximum height error, especially with MYJNoahNOURB (180 m; Fig. 7b, black line).Furthermore, using the MM5-TD LSM (Fig. 7b, orange line) allows for the correction of the day peak time advance.
During nighttime both parameterizations underestimate the height of the PBL and poorly reproduce the signal variability.The ability of both YSU and MYJ to capture the increase of turbulence during early hours of 7 February is important to note.The PBL growth during the morning of the second day is also overestimated by both PBLs (Fig. 7).The MYJ (Fig. 7b) anticipates PBL growth on 7 February by approximately six hours; nevertheless it attempts to reproduce the bimodal structure recorded by the lidar, even producing an overestimation of both the value and the duration of the maxima.It is anyway worth to note that mechanical contributions in MYJ simulations, even overestimated, act to suppress the thermal growth, as also found by Martilli (2002) during daytime.This is particularly true for the MYJNoahUCM2 simulation (Fig. 7b, pink line).On the other hand, the nonlocal-parameterization YSU (Fig. 7a) develops a typical diurnal growth of the PBL also during the second day, producing a large overestimation of the PBL height and allowing one to assess an excess of the thermal contribution to the layer growth.During the afternoon (13:00-18:00 UTC) both YSUtd (blue line) and YSUNoahNOURB (yellow line) rapidly decrease the layer height turning in its underestimation of about 300 m; in the following hours a new increase of the layer, clearly due to mechanical contributions, is produced, but again overestimating nighttime height measured by the lidar.Improvements are found by using UCM1 (Fig. 7a, green line) during 7 February, except for the evening overestimation.
The vertical structure of the PBL is further investigated using horizontal and vertical wind recorded by sodar (Figs. 8,9,10 and 11).The comparison between modeled and observed horizontal wind speed shows that WRF reproduces the dynamics occurring during the two days (Figs. 8 and 9) fairly well.The model overestimates wind intensities and anticipates the wind increase of the second day.The two PBLs produce similar correlations with observed data (0.84 for YSU; 0.83 for MYJ), but differences are found between them.The YSU produces a larger overestimation at lower levels than MYJ between 11:00 and 17:00 UTC on 7 February.The MYJ better agrees with sodar data for both the vertical profile variability during the first part of the simulation and for the development of wind speed maxima at higher altitudes than YSU (Fig. 8).This allows for weaker wind below 100 m than is represented in the nonlocal scheme and for a larger agreement with both sodar profiles (Fig. 8a) and observations inside the canopy layer (Fig. 5a).Moreover, an upward displacement of the maximum is produced by MYJ if UCM2 is activated (Fig. 9d), whereas no remarkable differences are found for YSUNoahUCM1 (Fig. 8d).The UCM2 activation correctly reduces the upper levels' air intrusion, decreasing the downward transmission of horizontal momentum.
The structure of the time series of the vertical profile of the model horizontal wind supports the hypothesis of a strong link between the canopy layer and the upper levels, causing the wind overestimation by the model with respect to the anemometer during second day, except for MYJNoahUCM2 (Figs. 3a and 5a,blue,yellow,green,orange and black lines).This simulation shows as an UCM acts to decouple the two layers (Fig. 9d), favoring a reduction of the forcing from upper to lower levels and thus producing a better agreement with the anemometer inside the urban canopy layer during the moderate advection regime (Fig. 5a, pink line).
The vertical wind components detected by sodar and WRF are shown in Figs. 10 and 11.To simplify the comparison, their color bars are not the same because of the large model underestimation.The model underestimates vertical motions regardless of the PBL parameterization, and poor correlation with measurements is found.In addition, WRF shows a longlasting downdraft, which is not detected by sodar, during the last part of the simulation.The YSU simulations do not show remarkable differences among them: the first day updrafts are missed by the urban canopy schemes, whereas only one updraft is produced by YSUtd at 12:00 UTC.During 7 February, upward motion is produced by all YSU at 06:00 UTC and between 10:00 UTC and 13:00 UTC.The differences between the MYJ simulations are found to be larger than those between the YSU simulations, and an overall better agreement between sodar and MYJ results is found than that between sodar and YSU results.Indeed, during 6 February the two updrafts are well reproduced using UCM2, with an even higher variability than the other MYJ simulations, but a delay with respect to sodar is found.During 7 February it improves the results by enhancing variability, but still anticipating maxima.
In summary, a different behavior in developing dynamics around the urban area of Rome is found using the two PBL parameterizations for this event.The results infer the different roles played by the local forcing during 6-7 February.Both days were characterized by clear sky, but, as already pointed out, during 6 February only weak wind was recorded, whereas during the following day, moderate wind advecting cool air was recorded.The model's ability in reproducing the first-day plume if using local PBL highlights the major role of the local forcing driving the dynamics; during the following day, both local and nonlocal forcing drive the dynamics; this induces difficulties for the model in reproducing the event in both configurations.Indeed, the model is not able to combine the two forcings; by balancing the large-scale forcing that tends to prevail, increasing errors occur with respect to observations.

The rural area
A further comparison between WRF and ground meteorological stations in the neighborhood of the urban area is performed.Nineteen of the twenty ground-based stations from SIARL agency 1 are used to investigate the effect of the different PBL on the local and regional circulation (the Lanciani 1 Servizio Integrato Agrometeorologico della Regione Lazio station located in the urban area has been excluded, Fig. 12).The SIARL stations recorded temperature (T2), relative humidity (RH), horizontal wind speed (WSP) and direction (WDR) at 2 m of height.High correlation values are found between WRF and observations for all meteorological parameters (Table 3) except wind direction.The YSU shows a slightly greater agreement with T2 and RH than MYJ, with YSUtd achieving the best values for all variables.The local PBL (MYJ) shows lower correlations than nonlocal one, with MYJtd having the highest correlation.Based on the correlation results, the YSU simulations show a better ability in reproducing the rural circulation than MYJ ones, whereas the opposite is found regarding the urban circulation.
The hourly averaged bias index, calculated as difference between observations and model results, shows a WRF tendency to underestimate/overestimate the observed T2 (Fig. 13a) maxima/minima, regardless of the PBL scheme and the land-surface model.Differences in the bias among the WRF simulations are found during the nigh time: YSU tends to underestimates the nocturnal temperature, whereas MYJ tend to produce unbiased results if the Noah LSM is used; the MYJtd still tends to overestimate observations in this interval (Fig. 13a, black, pink and orange lines).
Errors for RH (Fig. 13b) range between 0 and 30 %; a tendency of the model to underestimate maxima and overestimate minima has been detected.The MYJ shows a slightly better ability than YSU regarding the simulation of RH during nighttime.
The model shows a general overestimation of the horizontal wind speed (Fig. 13c), which increases after 24 h of simulation -that is, during the moderate advection regime.A poor sensitivity to all surface schemes is found for both PBLs except during the last part of the simulation.The overestimation with respect to the observations (0.23-0.98 m s −1 , using an UCM) are lower than the ones found in this study, but they are indicated as the lowest among other analogous studies.Flagg and Taylor (2011), for example, show a range of mean biases between −0.36 and 1.91 m s −1 using MYJ over Detroit with an accurate representation of the city area including different urban categories.The MYJ biases found in this work are in line with Flagg and Taylor (2011) (−1.49 to 1.24 m s −1 ) if we consider the MYJNoahUCM2 simulation (which is the most comparable in terms of configuration settings).
A cold bias is also often verified by other city studies.Kim et al. (2013) find a prevailing underestimation of temperature whatever the PBL and at different heights for Paris (mean value between −2.91 and 0.94 • C), with the bias decreasing with height.This underestimation is also confirmed by Wouters et al. (2013) over the same urban area (−0.92 • C).Underestimation of temperatures is also verified by Salamanca et al. (2011) and Flagg and Taylor (2011) for Houston (−1.46 • C) and Detroit (−4.90 to 0.38 • C) respectively.Underestimation found for the Rome area range between few tenth of degree to about 5 • C at the sonic anemometer height and to 2 • C at 2 m height in Lanciani station during nighttime.Sonic anemometer values are in the range of Flagg and Taylor (2011) findings in terms of maximum bias and the decrease of the temperature error with the height is in line with results in Paris (Kim et al., 2013).Finally, a study by Hu et al. (2010) also shows an underestimation of temperature, with MYJ producing larger bias than YSU, especially during nighttime, in agreement with findings of this work in the sonic anemometer comparison; Hu et al. (2010) address the underestimation to the surface layer scheme (Eta and Monin-Obukhov, MOY-MYJ and Monin-Obukhov, MO-MM5, respectively).
For the PBL height our results are in agreement with the findings of Kim et al. (2013) during daytime and in a thermaldriven regime: in both cases, a tendency to underestimate the layer height is found.On the other hand, opposite results are found during nighttime and under stable conditions for both MYJ and YSU.Further analysis is needed to justify this discrepancy.

Concluding remarks
In this study, an investigation of the circulation in the urban area of Rome is presented.Measurements from a sonic anemometer, a lidar, and a sodar, as well as the WRF model output, are used to highlight characteristics of the circulation in the urban area of Rome.An intercomparison of the two most commonly used PBL parameterizations (YSU and MYJ) is also performed with a twofold aim: evaluating the models' ability to correctly reproduce the PBL parameters, and understanding the different role of local and large-scale forcing in the area of interest.To investigate the impact of surface parameters on the PBL evolution, each PBL scheme is coupled with different combinations of a surface scheme together with a land-surface parameterization.
Results from numerical simulations at the highest resolution domain (780 m) are compared with measurements by the sonic anemometer, the lidar and the sodar inside the city.Moreover, a comparison between the model and observations for several rural-based stations is performed in order to investigate both the differences between the two PBLs and the impact of the urban forcing on the near-surface thermodynamical parameters in the rural area.The case study is selected based on the measurements availability and on its representativeness of a typical meteorological scenarios in the Rome area: 6-7 February 2008 is characterized by the transition from a weak large-scale regime to a moderate one, causing weak convection because of an increase of the wind speed.
Concerning the horizontal wind, the comparison with the anemometer reveals a tendency of the model to overestimate the horizontal wind intensity at low levels in both regimes, with larger errors if large-scale conditions prevail; the comparison with the sodar profile demonstrates the tendency of the model to develop wind maxima at lower levels than observed, thus suggesting an excess of vertical transport of horizontal momentum from upper to lower levels and an inefficiency in decoupling the canopy layer from the layer above.Only for MYJ using Noah land surface scheme and a multilayer urban canopy model, a decoupling of upper and lower layers is found.This results in a good reproduction of the horizontal wind field also in the canopy layer during the moderate wind regime, whereas a slight underestimation is produced during the weak conditions phase.
Regarding the vertical velocity, WRF largely underestimates sonic anemometer values regardless of the PBL parameterizations.Besides the large underestimation, a positive correlation with observed signal is found.
A tendency of the model to underestimate temperature is found for both regimes.YSU is generally able to correctly reproduce the evolution of the nocturnal cycle during the transition from weak to moderate large-scale regime.Time series show a poor sensitivity to different surface schemes, except for MYJ during daily cooling phases, when an increase of the error is found using Noah LSM.The relative humidity time series show a reduction of the error if MYJ is used especially during moderate wind regime.
PBL height retrieved from lidar shows a layer evolution mainly due to the thermal contribution during weak horizontal advection conditions, whereas mechanical contributions are also found during moderate/strong wind regimes due to large-scale circulation, which acts to reduce the PBL height with respect to the previous day.The comparison between the model and lidar revealed a tendency in the model to underestimate the PBL height if YSU is used for thermal prevailing conditions, whereas MYJ better reproduces this PBL structure.During the second day, when also a large mechanical contribution occurs, both PBL schemes overestimate the observations, with MYJ exceeding the mechanical contribution during diurnal hours but attempting to reproduce the signal variability.A poor sensitivity to the surface model is found, confirming that a major role is played by the mixing algorithms of the PBL parameterizations to the PBL height computation.
The comparison with rural area stations confirms the model's tendency to overestimate wind speed and produce larger errors if large-scale circulation is influencing the low levels.The temperature and humidity are generally underestimated/overestimated during daytime by both parameterizations.On the other hand, during nighttime YSU underestimates temperature whereas MYJ shows small biases.A poor sensitivity to the urban canopy model is found for both PBL parameterizations at the rural stations sites.
Based on this study, two major findings can be assessed: -the local forcing during the first day drives the local dynamics; whereas during the following day both local and nonlocal forcing drive the urban circulation, which is why the model, in both configurations, shows difficulties in reproducing the dynamics; -the MYJ scheme allows for reducing the errors of most of variables, even if further adjustments to the urban model parameters and to the surfaces schemes are necessary to improve results.
It will be of interest to achieve more definitive conclusions by statistically evaluating model performances over a longer time period in a future work.

Fig. 2 .
Fig. 2. Synoptic maps from ECMWF analyses at 0.25 • of resolution for (a) 6 February 2008 at 12:00 UTC and (b) 7 February 2008 at 18:00 UTC.Colors represent the mean sea level pressure (hPa), white lines the geopotential height at 500 hPa (m), and black vectors the horizontal speed at 10 m (m s −1 ).

Fig. 3 .
Fig. 3. Time series of horizontal wind speed (a), and direction (b), vertical wind (note: model result are multiplied by a factor of 10) (c), friction velocity (d) for 6-7 February 2008.Color code: YSUtd is blue; YSUNoahNOURB is yellow; YSUNoahUCM1 is green; anemometer measurements are red; and errors are grey.

Fig. 4 .
Fig. 4. Time series for temperature (a) and relative humidity (b) for 6-7 February 2008.Color code: YSUtd is blue; YSUNoahNOURB is yellow; YSUNoahUCM1 is green; measurements by probes combined with the sonic anemometer are red; and errors are grey.

Fig. 5 .
Fig. 5. Time series of horizontal wind speed (a), and direction (b), vertical wind (note: model results are multiplied by a factor of 10) (c), friction velocity (d) for 6-7 February 2008.Color code: MYJtd is orange; MYJNoahNOURB is black; MYJNoahUCM2 is pink; anemometer measurements are red; and errors are grey.

Fig. 6 .
Fig. 6.Time series for temperature (a) and relative humidity (b) for 6-7 February 2008.Color code: MYJtd is orange; MYJNoahNOURB is black; MYJNoahUCM2 is pink; measurements by probes combined with the sonic anemometer are red; and errors are grey.

Fig. 7 .
Fig. 7. PBL height time series for 6-7 February 2008 for YSU PBL (a) and MYJ PBL (b).Color code: YSUtd is blue; YSUNoahNOURB is yellow; YSUNoahUCM1 is green; MYJtd is orange; MYJNoahNOURB is black; MYJNoahUCM2 is pink; lidar-retrieved PBL height is red; errors are grey; and sodar-retrieved PBL height is a black, dashed line.

Fig. 8 .
Fig. 8. Time series of the horizontal wind speed vertical profile on 6-7 February 2008 for (a) sodar measurements, (b) YSUtd, (c) YSUNoah-NOURB, (d) YSUNoahUCM1.Time is on the x axis and height on the ordinate axis.

Fig. 11 .
Fig. 11.Time series of the vertical wind speed vertical profile on 6-7 February 2008 for (a) sodar measurements, (b) MYJtd, (c) MYJNoah-NOURB, (d) MYJNoahUCM2.Time is on the x axis and height on the ordinate axis.

Fig. 12 .
Fig. 12. SIARL stations (green points).The black circle approximately indicates the WRF urban area.The green point inside this area is Lanciani station.Red circle is the area with stations no more than 15 km far from the city.

. Pichelli et al.: Planetary boundary layer of the urban area of Rome as
the average between FDH and SDH (which are assumed to represent an estimation of the bottom and top heights of the entrainment zone, respectively).

Table 1 .
Outline of performed simulations.In the first column, the identification acronym of each simulation is shown.In the other columns, the parameterizations used are indicated (see the text for acronyms).
3 Meteorological analysis of the eventWith the aim of investigating the PBL structure in the urban area of Rome, a model-aided analysis of an event showing different meteorological scenario is carried out.An investigation of the model's capability to reproduce different meteorological scenarios is also presented.The 6-7 February 2008 event is chosen because local conditions associated with or without large-scale signals occurred in the urban area of Rome.To investigate the PBL structure over the urban area of Rome a comparison among the instruments is performed; the high-resolution model output will help in the understanding of local features.
). Color code: YSUtd is blue; YSUNoahNOURB is yellow; YSUNoahUCM1 is green; MYJtd is orange; MYJNoahNOURB is black; MYJNoahUCM2 is pink; lidar-retrieved PBL height is red; errors are grey; and sodar-retrieved PBL height is a black, dashed line.