The Micrometeorology of the Haifa Bay Area and Mount Carmel during the summer

: The Haifa bay area (HBA), which includes Mount Carmel and the Zevulun valley is the third largest metropolitan area in Israel. It is also a centre of heavy industry and an important trans ‐ portation hub which serve as sources of local anthropogenic pollution. Such sources are associated with adverse health effects. In order to estimate the possible exposure of the inhabitants in such heterogeneous orographic area, a detailed atmospheric transport and dispersion modelling study is required, which in turn must take into account the local micrometeorology. The aim of this study is to conduct a spatio ‐ temporal analysis of the flow field in the HBA in order to identify the common patterns of the average wind and characterize the statistical parameters of turbulence in this area, essential for detailed pollutants dispersion modelling. This study analyses data collected during four months of summer in a network of 16 weather stations which extend across Mount Carmel and the Zevulun valley. It was found that, during the evening and night time on Mount Carmel, different flow patterns may develop on each side, separated by the watershed line. When such conditions do not develop, as well as during the daytime, the wind field, both on Mount Carmel and the Zevulun valley is approximately homogenous. The analysis of the Monin–Obukhov similarity theory func ‐ tions for the velocity standard deviations show a distinct difference between Mount Carmel and the Zevulun valley, as well as between strong and weak winds. This difference can be clearly seen also in the diurnal hourly distribution of atmospheric stabilities which exhibit higher proportions of un ‐ stable conditions in the Zevulun valley during day time and higher proportion of stable stratifica ‐ tions at the Mount Carmel during night ‐ time.

the Mediterranean Sea, the eastern slopes of Mount Carmel and the Alonim range of hills [3].
Diversity also exists in the land cover in HBA. It consists of densely populated residential areas as well as major industrial facilities such as the national oil refiners, an oilfired power plant and several petrochemical, chemical and agrochemical industrial plants. In addition, several main roads pass through the HBA, especially the intercity roads that lead to Haifa, from the south (from Tel-Aviv), from the north-east (from the Krayot) and south-east (from Jezreel valley). These major industrial plants and the intercity roads serve as local anthropogenic sources of pollution, and especially particulate matter (PM) [4,5]. There are indications that increased PM concentrations may be associated with adverse health effects. Specifically, chronic exposure to ambient PM10 concentrations was found to be associated with lung cancer incidence among males in the Haifa bay area (HBA) [6]. In order to be able to examine the expected impact of different anthropogenic pollutant sources in the HBA on the health of people living in this area, a rigorous atmospheric transport and dispersion model should be used. Such a model should account for scenarios that correspond to the specific atmospheric conditions of this area [7]. Such scenarios should be based on comprehensive information on the average wind field, the statistical parameters of the turbulence as well as the influence of the heterogeneous canopy on these variables [8,9].
It was reported that the coastal plain of Israel is rather homogenous regarding dispersion of pollutants from power plants. However, this homogeneity excludes the HBA, due to presence of the complex topography of Mount Carmel [10]. Despite this fact, studies that dealt with the wind field in the HBA have based their analysis on very few weather stations, and did not include analysis of the time dependency of the wind [2,5,11].
Therefore, the aim of this study is to conduct a spatio-temporal analysis of the wind field in the HBA in order to identify and describe the common patterns of the average wind and turbulence in the HBA. It should be mentioned that the HBA is also affected by the regional semi-arid climate which is expressed by extreme natural dust events. As this study is supposed to provide information for models that will assess the effect of local emissions on the PM concentrations that people residing in this area are exposed to, this study will deal with the summer season, as during this season no natural dust outbreaks develop and anthropogenic sources dominate [12].

Weather Stations
The wind measurements analysed in this study originated from a combination of a dataset of a one-year wind measurement campaign conducted by the IIBR (Israel Institute for Biological Research) at the Haifa metropolitan area during 2014-2015. This dataset was supplemented by data from the ongoing monitoring network of Haifa metropolitan area which was installed and maintained by the HBUA (Haifa bay municipal association for environmental protection; http://www.envihaifa.net (accessed on 1 February 2021)). As this study deals with the summer season, the measurements that were made during the summer season were extracted from the dataset. As the summer season in the Eastern Mediterranean region spans from June to September, the data analysed in this study consist of the measurements collected during these four months in 2015 [13].
Overall, the weather stations' network analysed in this study consists of 16 stations, 10 in Mount Carmel and 6 in the Zevulun valley ( Figure 1). The wind measurements in the IIBR weather stations were equipped with ultrasonic anemometers (R.M. Young 81000), while the HBUA weather stations consisted of propeller-vane anemometers (R.M. Young 05103VM) ( Several field studies have compared co-located ultrasonic and mechanical anemometers, such as propeller-vane anemometers [14][15][16]. The main disadvantage of the later is over-speeding which occurs when an anemometer responds quicker to an increase in wind speed than to a decrease of the same magnitude [17,18]. An earlier study have estimated the over-speeding of a cup anemometer as 10% of the reference wind speed measured by a sonic anemometer, over a speed range of 1.5 5 m s ⁄ [16]. However, a more recent study that analysed a year of data have found that over a wider speed range of 0 12 m s ⁄ the differences were 2%. The over-speeding effect was noticeable in speeds up to 2 m s ⁄ , where the median bias was 0.13 m s ⁄ . Wind directions differences were within 5° 80% of the time [15]. Based on these results, the time averaged horizontal wind data from the weather stations network were used for the spatio-temporal analysis of the flow field in the HBA. In order to extract surface layer turbulence parameters only stations equipped with ultrasonic anemometer were used. Data quality control computer programs were applied on the data files, in order to remove erroneous records. The procedures included search for characters in numeric fields, examination of the time order of records, and examination whether the wind direction and speed were within the correct range, i.e., wind direction within 0°, 360° and wind speed within 0, 20 m s ⁄ . This maximal wind speed was set as the reported multiannual maximal wind speed in a representative station in Carmel coast was [11]. It should be noted that all reference to specific time is in LST (UTC + 2). All stations were installed on the rooftop of buildings. The height of the measurement above ground level (AGL) was between 5 and 30 m. The elevation above sea level (ASL) of the locations of the weather stations was between 0 and 350 m (Table 1).
Since the canopy roughness sublayer is highly inhomogeneous, in order for a station to represent its surroundings it should be just above the roughness sublayer, which is above at least 1.5 H (H being the averaged building height at the area of the station) for turbulent parameters and a little less for averaged wind if the canopy is sparse [19]. All stations in Table 1 matched the stricter criterion being in the range of 1.5 H-2 H, with the exception of station Z1, being at 1.3 H. This station is mounted on a 3 m mast on the rooftop of a single floor shelter (small cubed 1 floor structure of about 3m 3m 3m) where the surrounding a sparse urban area where the buildings are 1 (3 m) or 2 (6 m) floors.

Wind Field Interpolation
In order to analyse the spatial wind field, the wind measurements of the 16 weather stations were interpolated using a horizontal-and vertical-weighted wind field interpolation scheme that takes into account conditions of complex topography.
Two of the most popular interpolation methods in geoscience and atmospheric science is the inverse-distance weighting (IDW) and kriging, as indicated by the fact that they were implemented in many GIS packages [20,21]. Both of these methods estimate a value to a given point based on its distance to measurement points in space [21]. The distance that these methods rely upon does not take into account the vertical variability which is characteristic to complex terrains, and under a low sampling density may fail to describe important small scale phenomena, such as the difference between sunny and shady slopes or the difference between a valley and a nearby top of a hill [21,22].
In order to account for the effect of topography, and in particular where surface wind observations are concerned, it was proposed to replace the horizontal distance with the topographical elevation distance [22]. This study uses a surface wind field spatial interpolation method that was developed in the context of a detailed high resolution atmospheric transport and dispersion model [23]. This method follows the suggestion regarding the topographical elevation distance but instead of replacing the horizontal distance with the vertical, it incorporated both these dimensions. It is based on the IDW framework, due to the fact that the kriging is computationally cumbersome [20]. The method uses two length scales: a characteristic horizontal length scale, , ("influence" scale) by which the topography varies, and a characteristic vertical length scale, √ , associated with the influence of vertical atmospheric turbulent mixing. Therefore, the approach presented in this study expand previous methods for wind field interpolation by taking into account the extent of influence dictated by the topography as well as the effect of turbulent mixing [22,24]. The interpolation was performed on a square metric grid whose coordinates were in ITM (Israeli Traverse Mercator, a geographic coordinates system optimized for Israel). The wind vector values on a square grid of ∆ ∆ 100m were interpolated as a weighted mean of the values at the measurement sites (the weather stations' locations): where ⃗ ⃗ is the interpolated velocity at position ⃗ on the grid; ⃗ is the weight given to observation/measurement site for grid point ⃗; ⃗ is the observed velocity at site ; and the summation is performed over all measurement sites. The weight factors are taken as a product of two functions: First, a horizontal weight function , roughly based on the inverse square of ‖ ⃗ ⃗ ‖, the horizontal distance between grid point ⃗ and measurement site : where , here taken as 100m, is a characteristic horizontal length scale above which the topography starts varying considerably. Then, a vertical weight function , roughly based on the inverse square of now Δ , ⃗ ⃗ , the difference between the elevation at grid point ⃗ and measurement site : where is the same as above and is a dimensionless constant which depends on the atmospheric stability, taken here for the near neutral case to be 5 (for stable ∼ 100 and for convective atmosphere ∼ 2). The resulting characteristic vertical length scale over which information in transferred between different topographical elevations (by turbulent mixing) is accordingly /√ 45m.
This above interpolation is used in this study for identifying the typical diurnal wind patterns. The robustness of this interpolation method was tested by removing one station at a time (from the stations in Table 1) and comparing the interpolated grid calculated with and without this particular station. The error was defined as the influence of the omitted station to the spatial wind field, i.e., as the averaged L2 norm distance between the two grids, normalized to the average measured wind speed: where is the error obtained by removing station , ⃗ ⃗ is the interpolated velocity at grid point ⃗ obtained using all measurement stations, ⃗ ⃗ is the interpolated velocity at the same grid point obtained without station , n is the number of grid points and is the average over all stations of the measured wind speed.
The median error was found to be 3.4%, with an interquartile range (IQR): 3%-5%. We also note that the basis for choosing stations for sub-areal climatological description is that their error is being well within the IQR, as manifested in the choice of stations C4, C9, Z2, Z4 in Section 3.

Nonparametric Bivariate Density
The diurnal time series of the wind speed and direction in a given weather station may change from day to day. However, it is possible to identify climatological patterns in these time series, if each pair of time and wind vector component observation (either the direction or speed) is defined as a realization of a bivariate distribution. In this study, all these realizations were grouped into a two dimensional histograms, which were later plotted in a form of density contours. This resulted in a nonparametric visualization of the bivariate (time-speed or time-direction) distribution. Each contour represents a boundary in which a given percentile of the nonparametric bivariate distribution is contained.

Surface Layer Turbulence Parameters
Turbulent structure in the surface layer is usually described in terms of Monin-Obukhov similarity theory (MOST) using the well-known logarithmic velocity profile law and the so called universal relationships between scaling parameters [25]; The applicability of this description was critically reviewed [26]. The basic reasoning to MOST lies in two main observations. The first is that above the roughness sub-layer which extends to 2-5 averaged canopy elements height, the (horizontally homogeneous steady state) turbulent kinetic energy (TKE) budget approximately reduces to a "local equilibrium" between shear production and dissipation [27]. The second is that many times in this layer there exists a "constant-flux layer" in which fluxes are approximately constant with height (typically between 10-20% [28]). Analysis of the standard deviations of the wind-velocity fluctuations over complex surfaces is essential to our understanding of the turbulence structure of the atmospheric surface layer. In our case the terrain is characterized by varying topography and roughness. We focus here on characteristics of atmospheric turbulence over coastal flat terrain and a neighboring complex terrain under various wind directions, and atmospheric stabilities. These two cases allow us to compare the influence of terrain and roughness changes on the properties of turbulence, and in particular on the dimensionless velocity standard deviations, * ⁄ ( , , being velocity terms, and * being the friction velocity), which are central to the characterization of surface layer turbulence structure.
Turbulence data for the three wind components and the sonic temperature were measured at 30 Hz during the summer of 2015 (Table 1)

Results
The examination of two typical summer days shows that the night time in the Zevulun valley is characterized by low winds with easterly to south-easterly directions. During this time, two different regimes may develop in Mount Carmel (Figure 2a,b). One is characterized by south to south-south-west wind directions in the western slope of Mount Carmel and west to south-west winds in the eastern slopes of Mount Carmel (Figure 2a). The other is relatively uniform with south-westerly winds (Figure 2b). During the daytime, the wind directions measured on Mount Carmel are mostly west to west-southwest and in Zevelun valley west to north-west (Figure 2c,d). During the evening, the wind in Zevulun decreases in speed. The measured wind direction may exhibit a pattern similar to the daytime (Figure 2e), or a veering to the north-west (Figure 2f). During this time, the wind in Mount Carmel may be uniform, with south-westerly winds (Figure 2e) or a specific pattern in each slope, north-north-east directions in the western slope and northnorth-west in the eastern slope (Figure 2f). The wind direction in the Zevulun valley stations exhibit a similar behaviour during the daytime and evening (0700 h-2200 h). The wind direction distribution is uni-modal where its centre is around westerly winds. It is spread over a relatively narrow sector of 90°. However, while the central tendency of the wind direction distribution in station Z4 is 270° (West; Figure 4a), station Z2 experience some wind veering from 270° during 0700 h-1400 h to 315° (south-west) from 1400 h-2200 h (Figure 4b). The night time can be divided into two periods (2200 h-0300 h and 0300 h-0700 h). during the first night period, the wind direction distribution is multi-modal, where one mode continues around westerly directions (270°-315°) and the other is around 112.5° (East-South-East) in station Z4 (Figure 4a) and 22.5° (north-north-east) in station Z2 (Figure 4b). During the second night period the wind distribution in station Z4 is uni-modal around 112.5° while in station Z2 the wind direction is rather undefined, with modes around 112.5° (east-south-east), 157.5° (south-south-east) and 292.5° (west-north-west).
In the Mount Carmel stations, the central tendency of the wind direction distributions exhibit a clockwise veering from the south-east to the north-west (Figure 4c,d). This occurs during the daytime (0800 h-2000 h). During most of this period (until 1600 h) the wind direction distribution is uni-modal. From then on, the distribution is multi-modal, with modes in the west, north-west (270°, 315°, 337.5°).
During the night time, there is a distinct difference between the Mount Carmel stations. Station C4, which is located in western slopes of Mount Carmel, exhibit two separate periods. During the first (2000 h-0400 h) the wind direction distribution is multi-modal, with modes in the North-West, West and South-West (315°, 270°, 225°). The second period is characterized by a uni-modal south centred distribution (0400 h-0800 h; Figure 4c). Contrary to this wind regime, the night time wind distribution in station C9, on the eastern slopes of Mount Carmel is uni-modal during the whole period (2000 h-0800 h). First, the distribution is western centred (270°; 2000 h-0400 h) and then south-eastern centred (225°; 0400 h-0800 h; Figure 4d).
Compared to the wind direction regimes in Mount Carmel and Zevulun valley, the wind direction regime in the Tel-Aviv area exhibits four distinct periods. During the day, from 0800 h, winds that rotate clockwise from westerly to north-westerly winds. During the night, from 2200 h in TO and 2330 in TI, south-westerly winds. Between these periods there are transition periods (Figure 4e,f). During the transition periods the wind changes from the night to day regime and vise-versa. The rotation of the wind may occur either clockwise or anticlockwise [29].
The analysis of turbulence parameters was based on the data measured at station Z2 (Kiryat Motzkin) which represented the flat urban region that is the Zevulun valley, and C9 (Tel Hai elementary school) that represented the complex terrain ridge of Mount Carmel (Figure 1 and Table 1). The measurement sensors at both stations were stationed at about twice the mean height of canopy elements ( ), which for closely spaced (urban) canopies is assumed to be above the roughness sub-layer [30,31].
In Figures 5-7 the dimensionless turbulence parameters for the vertical and horizontal standard deviations of wind-velocity fluctuations / * are plotted as functions of z/L. Similar to [32] we have evaluated turbulence parameters on the available data for the full range where * is the friction velocity calculated in all cases from formula * , is the stability parameter, z is the height, * is the Obukhov length (g being the gravitational acceleration, and 0.4 the von Karman constant), and d is the zero-plane displacement height, taken here as 0.7 [33]. As can be seen from the figures, the observations around the neutral stratification, within the range of 0.1, show constant value, after which it grows on both sides. The fit to the empirical coefficients , is given in Table 2.  In Figure 8 the diurnal hourly distribution of stabilities for the summer season based on the stability parameter is presented. The five categories we chose are adopted from [34]: for the neutral |1⁄ | 0.0081, very stable/unstable stratification | 1⁄ | 0.875, and in between the stable/unstable regime. A very distinct difference can be seen as the station on Mount Carmel exhibits a more pronounced stable stratification and less convective stability compared to the station in the Zevulun valley.

Discussion and Summary
This study has analysed the diurnal and spatial patterns of the wind field in the HBA, with its two distinct regions of Mount Carmel and Zevulun valley. It was found that the examined weather stations in this region exhibit similar wind speed diurnal regime. In this regime the wind speed is weak (~0.5-2 m/s) during the night and grows stronger during the morning until a maximum of around 3-6 m/s is reached around noon time. It was also found that this regime is similar to metropolitan area of Tel-Aviv. This result is in line with the wind speed regime farther to the south of Israel, the southern coastal plain [35]. As the HBA, metropolitan Tel-Aviv and the southern coastal plain are all regions that are adjacent to the Mediterranean Sea, it can be deduced that they are all affected by the landsea temperature difference. This phenomenon is closely related to the sea breeze [35].
The spatial analysis of the wind field in this study has shown two different flow regimes that may develop during the night time in Mount Carmel. One regime is characterized by a relatively uniform south-westerly winds. The other manifests a divergence of the flow between the two slopes, where on the western slope south to south-south-west wind directions occur and on the eastern slope west to south-west winds. This phenomenon of two wind regimes, a uniform one and the other separated by the watershed line was also found in relation to wintery rain events in Mount Carmel. The uniform regime was characterized by southeastern winds, while during the separated regime, mainly south-western to western winds were measured on the western side of the mountain, and on the eastern side northwesterly or southeasterly winds [2]. The fact that this and the previously mentioned studies, although each deals with different seasons, and finds different typical winds, find flow patterns of either uniform or separated flow patterns on Mount Carmel, supports the insight that the topographical disturbance posed by Mount Carmel may occasionally lead to different flow patterns in each of the mountains' sides.
The analysis of the diurnal regime of the wind directions in the HBA has shown periods of uni-modal wind direction distribution during the night time, although the wind speed is rather weak. During these periods, the wind distribution is centred around the south or south-east on Mount Carmel and mostly eastern or south-easterly in the Zevulun valley. These periods begin on 0400 h on Mount Carmel and 0300 h in the Zevulun valley. Compared to these, south-westerly winds develop in outskirts of Tel-Aviv from 2200 h and in inner Tel-Aviv from 2330 h. In the case of Tel-Aviv area, it is possible that this area is influenced by down slope flows from the Judean mountain range.
The night time periods of uni-modal wind directions' distributions may be associated by the land breeze that develop due to katabatic winds and land-sea differential cooling. Interestingly, a numerical weather prediction (NWP) model with a space-time resolution of 2 km-6 h that was used to reconstruct the timing of the land breeze phenomenon in the Israeli summer was able to provide a rather accurate timing for the Tel-Aviv area (2100-0100 h), but predicted a too early beginning of this period (2300 h) in the HBA [36]. This suggests that the complex and heterogeneous topography of the HBA requires a specific modelling effort in order to be described accurately. Figures 5-7 show a fit to the semi-empirical MOST functions (Equation (5)) for the stable as well as the unstable regimes. We should note that the free convective regime is known to scale using mixed layer height [37]. According to [38], below the free convective layer there exists two more sublayers, a shallow dynamic sublayer up to ∼ 0.3 and a convective-dynamic sublayer above up to ∼ 3 . Both our measurements are at about twice the averaged canopy height, and therefor are expected to be mostly in the convective-dynamic sublayer (which scales with ).
The analysis of turbulent parameters in complex terrain, such as the respective standard deviations of the longitudinal, transverse and vertical components of the wind, might be affected, especially for developed winds, by an upwind change of surface roughness (e.g., [39,40]). This issue can be addressed by analysis according to wind direction sectors correlated with different upwind roughness (e.g., [39]). The wind characteristics in the HBA, as exhibited in Figure 4, are essentially confined in a single sector during night (with weak wind and stable to neutral stratification) and a single one during the day (with developed winds and unstable to neutral stratification). Therefore, there is no need to address sectors of difference surface roughness.
Regarding the values of our empirical fits (Table 2) to the parameters of Equation (6), these are well within the common range of values for these parameters documented the literature (e.g., [26,41]), although in the ridge station horizontal homogeneity is only approximate. We should note that in accordance with [39] we also find that there is a distinct difference between data points measured in the valley and the ridge stations. In the neutral stability limit ( / → 0) we obtain 1.25, 1.15 for valley and complex terrain, respectively. This value agrees with the mean value of given by [42], which is 1.25 over flat terrain and is lower than their 1.23 over complex terrain, and in accordance with the value of 1.2 of [39] for a station surrounded by different topography. Our valley values are lower than those of [43] which are around 1.3 and 1.4 for along-and cross-valley, respectively. Our interpretation to that is that the Zevulun valley's station in this study resembles a flatter area than a typical valley.
In contrast to the values of , the first parameter of Equation (5), the estimates for the second parameter, , vary dramatically between strong ( 1 m s ⁄ ) and weak ( 1 m s ⁄ ) wind speeds and between the two sites ( Table 2). The difference between strong and weak winds is especially pronounced in the valley station in convective conditions; due to the urban heat island phenomena however, the valley station recorded too little data in the stable regime for the difference between strong and weak winds to be statistically significant-and this is further corroborated by the low frequency of stable regimes in Figure 8 discussed below. In the ridge station, a difference between strong and weak winds was only found for the vertical standard deviations. From the hourly distribution of stabilities for the summer season (Figure 8), both stations, the flat urban valley station and the complex terrain ridge station, display neutral and stable conditions during night-time, with neutral ones being more common. In the flat valley urban station, we also have a small contribution of unstable conditions. During the day both unstable (convective) and neutral condition are dominant, with the fraction of neutral conditions significantly increasing during the day, becoming pronounced from early afternoon till evening. However, it can be clearly seen that these two stations, exhibit pronounced differences. In the ridge station the unstable and certainly the very unstable day conditions are less pronounced. This can be attributed to the stronger winds typical to mountain top ridge areas. Regarding the night-time, the ridge station shows a higher proportion of stable stratifications, which can be attributed to two contributions: the heat island forming more effectively in flat areas, and the higher elevation or shallower atmospheric column for the ridge station leading to more efficient radiative cooling.