1Introduction

Nowadays, it is well established that the thermal internal boundary layer (TIBL) is formed by the sea breeze. Many studies have noted and investigated the presence of this layer (Garratt and Physick, 1985; Physick et al., 1989; Challa et al., 2009). From a practical point of view, the TIBL traps pollutants released in the lower altitudes via the fumigation effect and can reintroduce pollutants present above, increasing the pollutant concentration in areas near large volumes of water (Bouchlaghem et al., 2007; Boyouk et al., 2011). Therefore, a realistic reproduction of the TIBL by a mesoscale model is crucial to improving the accuracy of subsequent air quality modeling and assessment.

Sea breeze is a well-known mesoscale meteorological phenomenon (Stull, 1988). Due to the daytime heating of the ground by the sun, the portion of air just above the surface is heated, expands, and rises. Consequently, a low-pressure region is formed. Combined with light surface winds, this area of low pressure can be occupied by cooler air from the sea surface to form a circulation cell directed towards the land at lower altitudes and towards the sea at higher altitudes. At night, if the temperature over the land becomes less than that over the sea, then the process is reversed to form a land breeze. Though seemingly simple, the intensity of these phenomena is very complex (Melas et al., 1998; Miller et al., 2003; Porson et al., 2007) and depends on several variables such as soil moisture, which consumes energy in the evaporation process, topography, soil type, land use, and speed and direction of the synoptic winds.

Surface characteristics determined with the use of land surface models (LSMs) are vital because the incidence of sunlight at the surface initiates the physical processes that subsequently define the structure of the atmospheric boundary layer (ABL) (Xiu and Pleim, 2001). Such processes cannot be described by the traditional equations of fluid mechanics. Therefore, parameterizations that can effectively describe the spatial and temporal evolution of these processes have been studied by several researchers (Chen et al., 1996; Koren et al., 1999; Chen and Dudhia, 2001). Experimental studies and numerical approaches have also been used to understand the ABL structure in regions that experience sea-breeze phenomena (Finkele, 1998; Miao et al., 2003). For example, in the work of Papanastasiou et al. (2010), the Weather Research and Forecasting (WRF) model was employed to study the wind field over the east coast of central Greece under typical summer conditions. The results from the model simulations were compared to observations that were collected by a near-surface meteorological station. In the work of Muppa et al. (2012), the characteristics of sea-breeze circulations over a tropical Indian station were studied, based on one year of observations, by using SODAR (sonic detection and ranging). Other related studies include Kitada (1987), Clappier et al.(2000), Talbot et al.(2007), and De LÉon and Orfila (2013). Several studies have devoted to the evaluation of WRF parameterizations in sites not subject to sea-breeze phenomena (Zhang and Zheng, 2007; Borge et al., 2008; Han et al., 2008; Li and Pu, 2008; Hu et al., 2010; Ma et al., 2011; Shin and Hong, 2011; Cheng et al., 2012; Osuri et al., 2012; Xie et al., 2012; Garcia-Diéz et al., 2013; Hu et al., 2013). However, the literature lacks contributions that evaluate different ABL and LSMs schemes in seabreeze conditions with simultaneous LIDAR (light detection and ranging) and SODAR high-frequency observational data.

In this study, experimental and numerical approaches are combined to study a sea-breeze system in a flat coastal area of the North Sea. The main objective is to study the predictive abilities of parameterizations used by the WRF model (version 3.4.1) to simulate a sea-breeze system in Dunkerque, France, with special attention to determine the TIBL height due to the problems with air pollution in this region. These include Mellor-Yamada-Janjic (MYJ) (Janjić, 1990, 1994), Yonsei University (YSU) (Hong et al., 2006), and Asymmetric Convective Model version 2 (ACM2) (Pleim, 2007a, b) parameterizations for the ABL; and Noah (Chen and Dudhia, 2001) and the rapid update cycle (RUC) (Smirnova et al., 1997, 2000) LSMs. The prognostics of the sea-breeze dynamics output from the WRF meteorological model were compared with measurements taken by SODAR, LIDAR, and MSS (meteorological surface station) to verify that the model accurately predicted the behavior of local circulations.

The paper is organized as follows. In Section 2, the details of the measurements are described and the model description and initialization used in the simulations are introduced. In Section 3, the meteorological monitorings by using MSS, LIDAR, and SODAR are presented, including comparison with modeling of the sea breeze. Then, a brief discussion is presented in Section 4. Finally, in Section 5, the conclusions are provided.

2Experimental set-up and model 2.1Measurement campaign

The experimental campaign was conducted in a French industrial coastal region on the North Sea (see Fig. 1a), in Dunkerque, from 7 to 9 September 2009. However, the sea breeze was present only on 8 September. Heavy industrial activities in this coastal area in the northern part of Dunkerque generate hazardous pollution episodes in the urban zones, which are subject to significant industrial and maritime influences, especially during a sea-breeze circulation (Talbot et al., 2007; Rimetz-Planchon et al., 2008).

Six meteorological parameters measured from the SODAR, LIDAR, and MSS (the MSS is located at atmospheric monitoring organization-the red point in Fig. 1e) were used to evaluate the WRF model, i.e., temperature, heat flux at the surface, horizontal and vertical wind speeds, wind direction, and ABL height.

The observations were made by using groundbased remote sensing instruments and a meteorological surface station in the Dunkerque area (Fig. 1e). To study the dynamics of atmospheric variables, a UV scanning LIDAR (Fig. 1b), a wind profiling SODAR (Fig. 1c), and an MSS (Fig. 1d) were deployed near the industrialized coastal area. During the campaign, an angular aerosol LIDAR (Leosphere ALS 300) was operated with the third harmonic of an Nd-YAG laser (355 nm), which has an energy pulse of approximately 16 mJ at a repetition rate of 20 Hz. High spatialtemporal resolution vertical and horizontal scans with continuous monitoring were conducted to obtain the vertical structure at lower altitudes and the dispersion of pollutant smoke plumes over a range of approximately 5 km. The LIDAR was installed on the roof (12 m) of the MREI (Maison de la Recherche en Environnement Industriel) building at approximately 1 km from the coastline. The LIDAR’s blind distance for optimal near-field overlap was 250 m, and the spatial resolution was 15 m along each beam. The acoustic Doppler SODAR system (Remtech PA2) was located at MG (Marégraphe) along the coastline and close to a 15-m mast equipped with an ultrasonic anemometer. Despite its specified maximum height of 1000 m, it provides valid information on the vertical wind profile (wind speed, turbulence, vertical movements) from 25 to 500 m (vertical resolution 25 m, time interval 15 min). Data for speed and wind direction were obtained at 10 m, and temperature was measured at 2 m by the MSS (ultrasonic anemometer) (Fig. 1d).

2.2Synoptic conditions

During the first two days of the experimental campaign, a high pressure system presented over North France, dominating the weather with stable conditions. Anticyclone air descended, forming an area of high pressure at the surface. Because of these stable conditions, the weather was settled, with only small amounts of cloud cover. These synoptic conditions were conducive to sea-breeze development during the field campaign. At the end of the campaign, a cold front was over the Dunkerque region, making the weather unstable.

Calm and moderate winds were measured on September 8 and 9 at 0000 UTC (Local Time is UTC + 1), showing a maximum value of horizontal wind speed at approximately 1500 UTC (6 m s^-1). When the cold front was over the Dunkerque area (September 9 to 10), strong winds occurred, and wind speed increased above 10 m s^-1. During the pre-cold-front time, on September 8, the wind comes from SE-S-SW until noon. After that, the main wind comes from W-NW until night, showing the effect of the sea breeze in this region. Due to the influence of the cold front, on September 9, the wind was almost constant during the day and night, and comes from N-NW. Dunkerque had high temperatures throughout the field campaign, varying from 16℃ (night) to 26℃ (day). During the cold front passage (cooler air), Dunkerque had cloudy skies and an average of temperature approximately 18℃; it showed little variation from day to night.

2.3WRF model set-up

As shown in Fig. 2, the WRF model was run from 6 to 10 September by using four nested domains (D1, D2, D3, and D4). D1, D2, and D3 are squares, with sides measuring 1863, 891, and 297 km, respectively, and 27-, 9-, and 3-km grid cell spacing, respectively. The inner domain (D4) represents the area of our study, with sides measuring 63 km in the W-E direction, 81 km in the S-N direction, and 1-km grid cell resolution. The four domains were centered at Dunkerque City, France, at 51.035°N, 2.376°E (the same coordinates where the measurements were obtained) and had 43 vertical levels, up to 13000 m, with more levels concentrated near the ground (32 levels under 1600 m). The height of the geopotential center of the lower cell started at 2 m and had a spacing of 25 m until the 500-m altitude. Above this point, the grid spacing increased because there was no SODAR data.

To evaluate the results predicted by the WRF model simulations, observations were made with several instruments, including SODAR, LIDAR, and MSS. The six combinations of physical parameterizations (ABL and LSMs) that were used in the WRF model are shown in Table 1.

The simulation was performed by using an oneway mode with a 54-s time step in an external domain and a reduction ratio of 3 for the other domains. To set up the initial and boundary conditions, NCEP FNL Operational Global Analysis data were obtained, with a 1° spatial resolution and a 6-h temporal resolution. Terrain/land data were supplied by the U.S. Geological Survey with 10-, 5-, and 2-min resolutions for D1, D2, and D3, respectively, and 30-s resolution for the innermost domain. Other parameterization (microphysics and radiation) remained unchanged for all of the runs (WRF single-moment 6-class scheme to microphysics; rapid radiative transfer model for longwave radiation, and Dudhia scheme for shortwave radiation). Due to the grid dimensions, the parameterization of cumulus was an exception, with the use of Betts-Miller-Janjic for the first and second domains, Grell-3D for the third domain and off in the fourth domain because it had a grade lower than 3 km. A detailed description of the basic equations and numerical schemes of the WRF model can be found in Skamarock et al. (2008).

3Numerical results

In this section, to verify that the model has reasonable accuracy in predicting the sea-breeze dynamics, the WRF model simulations were compared with the measurements from MSS, LIDAR, and SODAR.

A LIDAR system measures light scattered by the atmosphere, so the backscattered signal contains information about the extinction and scattering of light by aerosols. The aerosol load differs for each atmospheric layer and front; thus, the derivative of the LIDAR signal, which is called the negative LIDAR signal variation (NLSV), is discontinuous between the layers and provides information about atmospheric structure (see Melfi et al., 1985). A boundary between layers could be detected by the inflexion point method (Menut et al., 1999), which finds the position of the minimum of the second derivative of the backscattered power measured by the LIDAR.

In Fig. 3, the NLSV LIDAR data (in color) are coupled with the inflexion point method employed for ABL and TIBL detection (black diamonds). This image shows the ABL growth and the sea-breeze emergence beginning at approximately 1300 UTC to form the TIBL.

Figure 3 shows the beginning of the ABL growth at approximately 0700 UTC. At approximately 1300 UTC, it droped to less than 250 m, and there was a significant increase in the LIDAR signal intensity (a.u., arbitrary units), which was associated with the growth of the TIBL. It is noteworthy that, between 1200 and 1300 UTC, the ABL had an upward trend. However, the LIDAR detected the presence of this new region below 250 m, associated with the emergence of TIBL, which caused a decrease in signal intensity in the upper regions. Between 1600 and 2000 UTC, at altitudes between 250 and 1000 m, the signal maintained a greater intensity, and at approximately 1600 m, the signal became very intense, characterizing the presence of a higher region associated with the ABL. At approximately 1600 and 1700 UTC, there was again a decrease in the signal below 250 m, and it increased again at approximately 1800 UTC until 1930 UTC (suggesting a TIBL with lower intensity). At that point, the signal again decreased, suggesting the emergence of a new region below 500 m from 2000 to 2400 UTC.

The sea breeze depends on a large set of factors that can substantially alter its behavior, both in spatial and temporal terms. As the soil near the shore is heated by sunlight, the local air parcel rises and expands, creating a low pressure zone. Under the right conditions, this zone can be occupied by the displacement of a portion of the air that is stable over the ocean. This creates at low altitude, colder air current (stable) unlike that of the heated air (unstable) due to larger-scale phenomena. However, the cold air layer moving from the ocean to the land is also heated by the ground, forming another unstable layer near the ground (TIBL).

The ABL height is a crucial parameter because it represents the limited region of boundary layer heating, convection, and low-level circulation. The rapid decrease in the mixing depth (1300 UTC) during this measurement campaign resulted from the development of a TIBL associated with the sea breeze. A similar appearance of a new layer due to TIBL was observed during the Expérience sur Site pour COntraindre les Modèles de Pollution atmosphérique et de Transport d’Emissions experiment over the Marseille area in southern France in 2001 (Delbarre et al., 2005; Puygrenier et al., 2005).

3.1WRF and MSS data

Figure 4 shows a comparison of the MSS measurements and WRF horizontal wind speed (m s^-1) simulations near the ground for all runs.

All sets can reproduce the increase in synoptic wind speed until noon, with a maximum speed of the MSS data of approximately 5.3 m s^-1. Around the sea-breeze frontal passage (before 1300 UTC), the wind speeds from sets 1 and 5 had better agreement with the MSS data. However, sets 2 and 6 could also capture the wind speed drop at this time, with set 2 having a higher peak speed. Sets 3 and 4 could not reproduce this fall accurately. Regarding the second increase in wind speed after 1600 UTC, only set 4 agreed with the MSS data. Specifically, the synoptic wind speed directly influenced the further reduction of the breeze intensity.

Figure 5 shows a comparison of the MSS data and the WRF horizontal wind direction simulations near the ground for all runs. The shift in the wind direction can be used as an indication that a sea-breeze front passed through the MSS in the Dunkerque area. This behavior was observed during the experimental campaign when the predominant wind direction changed from S-SW to W-NW at approximately 1300 UTC.

Figure 4 shows an increase in wind speed (2.5 to 5.3 m s^-1) between 0800 and 1200 UTC. Figure 5 shows a shift in the wind direction between 1230 and 1400 UTC from southerly to northerly under the influence of the sea-breeze onset. The plots of sets 3 and 4 indicate that the model values differed from the observations due to a delay in the simulated wind direction. All runs failed to replicate the change in wind direction at 1700 UTC, i.e., the north wind directed S-W and returned northward in only one hour. Sets 1, 2, and 6 have good agreement with the wind direction at the sea-breeze passage time.

Figure 6 shows the temperatures obtained by the MSS compared with the model output temperatures for the grid cell containing the measuring point.

At night, the temperature at 2 m was significantly underestimated, and the humidity was always overestimated (figure omitted). Sets 3 and 4 showed delays relative to the measured values at the time of the seabreeze onset.

During the morning, set 4 showed the greatest temperature difference (6.9℃). All sets underestimated the temperature in the daytime, which was attributed to a persistent underestimation of the nighttime temperature. Sets 3 and 4 presented reasonable delays relative to the measured values at the time of the sea-breeze onset. Sets 1 and 2 represented the decrease in the temperature due to the breeze well, although sets 5 and 6 also had coincident times of temperature drop, however, they underestimated the temperature by almost two degrees at the moment of the fall. No scheme captured the temperature rise at 1700 UTC; the temperature instead decreased for all sets.

Analyzing the previous figures makes it possible to conclude that the second temperature peak at approximately 1700 UTC was associated with the reversal of the wind direction because, at that moment, the wind from the sea was slower than that from the land (synoptic wind) which carried warm air from the continent; this phenomenon was not simulated in any of the six sets. The breeze penetrated the continent, but the synoptic wind was in opposition. This behavior explains the results of the TIBL gaps from LIDAR and SODAR data in the sequence.

In numerical models, the air temperature at the roughness height (z₀) is often replaced by the surface radiative (skin) temperature, which can be computed from the surface energy budget. The radiative skin temperature calculated in the LSM, T=T_s, is used as a lower boundary condition (at z=z₀) for the surface layer parameterization, in which surface exchange coefficient for heat and moisture is calculated. Figure 7 shows the behavior of the skin temperature for all sets.

In set 6, the skin temperature was near 30℃ at approximately 1300 UTC. In this period of the inwarddirected breeze, at approximately 1400 UTC, the temperature decreased by 4℃, indicating sharp cooling of the surface. The temperature continued to decrease until nearly 1800 UTC and reached approximately 14℃. Compared to set 2, which changed only the parameterization of the surface, set 6 behaved similarly at 0600 and 1200 UTC and suffered a time delay in the temperature drop. The observed delays were correlated with skin temperature because, with the sea-breeze entry, these temperatures were maintained with little variation. They only decreased after a certain amount of time. Thus, the influence of the LSM scheme, which calculates the skin temperature, was clearly evident. Sets 3 and 4 had sharper drops in skin temperature that almost coincided with the temperature drop delay at 2 m because of the sea-breeze input.

Skin temperature is needed in calculating sensible and latent heat fluxes; specifically, sensible heat flux is determined by the instantaneous difference between the near-surface air temperature and land surface temperature (LST). Besides, LST also determines the amount of thermal heat that is vertically transported into the ground. As an attribute of the land surface, LST is influenced by the local land-cover.

The behavior of the sensible heat flux is shown in Fig. 8. During the day, the heat flux was overestimated for all sets, except around the sea-breeze start time. The role of LSMs is to redistribute the incoming solar radiation flux between latent and sensible heat fluxes, depending on precipitation and atmospheric surface layer parameters. According to analyses of sensible heat flux, RUC and Noah schemes were found to overestimate the sensible heat flux. This can explain the biases of these LSMs in our simulations. This might be due to surface parameters, such as land use, leading to problems of the model representativeness.

A significantly higher sensible heat flux was found in the experimental dataset at approximately 1300 UTC. In addition, a negative heat flux was observed during 1600-1800 UTC. In the morning (0900-1200 UTC) and afternoon (1400-1800 UTC), the model overestimated the sensible heat flux in all runs; the model did not predict the abrupt sensible heat flux changes, including the negative flux between 1600 and 1730 UTC. In particular, Sets 2 and 6 gave better representations of the experimental sensible heat flux peak (225 W m^-2) at approximately 1300 UTC. There is a relationship between the temperature at approximately 1700 UTC and the negative sensible heat flux from Fig. 8; at this time, this behavior is probably because the ground was colder than the air above. This low experimental sensible heat flux reduced the convection and turbulence, delaying the onset of the seabreeze circulation.

As stated previously, it is possible that the negative heat flux is also associated with the reversal of the wind direction. In fact, this suggests that the synoptic wind carried warm air from the continent due to the low intensity of the sea breeze, and the region above the ground became more heated. The soil was colder due to previous air inlet, so there was a negative heat flow at this time. However, this behavior was not observed in any of the six sets. The LSM information used in this study failed to capture these changes. It is relevant also to consider that the type of soil is very important, and this information probably was not efficiently introduced into the model with the horizontal grid of 1 km.

3.2WRF and LIDAR data

Compared with measurements from LIDAR, the most accurate estimation of the ABL was obtained by using set 6 (Fig. 9). In this run, the ABL height increased until 1245 UTC to 1000 m before the sea breeze began. The WRF model detected the mixed layer until 1000 m, when it started to detect the TIBL. The TIBL top could be monitored and was located between 100 and 200 m, from 1400 to 1545 UTC. After this time, the TIBL height began to decrease until late afternoon (Fig. 9).

The ABL height, particularly the TIBL height, is a key variable for air quality modeling. Its accurate simulation is often difficult in numerical models. Figure 9 shows that sets 1, 2, 3, 4, and 5 produced the largest ABL heights compared with ABL height from LIDAR, and sets 3 and 4 produced the greatest delay in detecting the sea-breeze passage. The local ABL scheme used in set 6 best predicted the ABL height, while the other sets overestimated the height by approximately 20% compared with the LIDAR prediction. Notably, the ABL height diagnostic formulation is specific to the ABL scheme. The formula for diagnosing ABL height (h) is specific for each parameter. YSU and AMC2 determine the top of the ABL where the bulk Richardson number exceeds a critical value, and the YSU design calculates the bulk Richardson number based on surface soil. Local designs of the ABL, of the MYJ type, define this height as where the turbulent kinetic energy decreases to a prescribed value. According to the results, the MYJ design shows a more satisfactory performance in determining ABL height in sea-breeze conditions in this region for the study day. Furthermore, the results suggest that the Noah parameterization could best capture the energy balance information in the soil.

Because the ABL height is determined by using different methods, it is better to make comparisons between the LSMs by using the same ABL parametrization.

Regarding sets 1 and 5 (YSU), there was similar behavior, with a maximum ABL height of approximately 1200 m. However, set 1 fell more quickly and with less delay, suggesting that RUC parameterization best captured the surface information.

Analyzing sets 2 and 6 (MYJ) shows that there was an underestimation of ABL height and an overestimation of the maximum height of the ABL (1400 m) at approximately 1300 UTC for the results of set 2, but it shows a delay for the TIBL. However, set 6 results showed a good agreement with the growth of the ABL, with a maximum of approximately 1000 m at 1130 UTC, ranging 800-1000 m from 1230 to approximately 1300 UTC, when the presence of the TIBL was detected. After, it shows a slight underestimation of the TIBL. Comparing sets 2 and 6, the results suggest that Noah parameterization gave better performance. Note that YSU and MYJ have the same surface layer scheme.

Regarding sets 3 and 4 (ACM2), there was fairly similar behavior, with a height of 1200 m, but these sets have the longest delays in time simulations of the TIBL detection. The longest delay is given by set 4, suggesting the poor performance of the Noah parameterization.

Thus, MYJ performed satisfactorily in determining the ABL in these sea-breeze conditions, and Noah LSM somehow managed to obtain better information from the energy balance in the soil. Both schemes, YSU and ACM2, did not satisfactorily capture the information from soil to avoid this delay. Figure 10 shows additional information about the TIBL and ABL obtained from potential temperature simulations.

Set 6 showed, at approximately 100-m altitude, a potential temperature gradient that is well defined at 1300 UTC, which did not occur for the other sets. For larger heights, above 1000 m, all sets have a potential temperature gradient to determine the ABL height. Clearly, it is possible to identify the coexistence of two inversions, or two regions at the same instant in time. The same observation was made to 1500 UTC, but sets 1 and 2 also presented a reversal at low heights.

3.3WRF and SODAR data

Figure 11 shows temporal comparisons of the horizontal wind speeds and wind direction. In the SODAR data, the surface winds turned from S-N at 1200 UTC to N-W at 1300 UTC, indicating the onset of the sea breeze (speeds of 4-6 m s^-1). However, all sets overestimated the observed wind speeds at high altitudes. This difference could be related to the observed daily evolution of the surface temperature (Dai and Deser, 1999). Temporal comparisons of the horizontal wind speed and wind direction reveal that all sets did not do well in detecting the sea-breeze frontal passage. Set 6 had similar wind speeds, but all other sets showed a delay relative to the frontal passage measured by the SODAR (change of wind direction at altitude). The delays between the measurements and the WRF outputs were as follows: set 1: 0.5 h; set 2: 0.5 h; set 3: 2.75 h; set 4: 3.25 h; and set 5: 1.25 h. Most of the runs represented the wind direction and speed fairly well at high altitudes; only set 6 of the WRF simulations showed a similar intensity and direction of wind as the SODAR observations on the ground level before and after the sea-breeze passage.

A better view of the SODAR data is shown in Figs. 12 and 13, respectively, comparing with wind speed and direction data. These WRF model simulating are for all sets at three different times: 1230, 1300, and 1330 UTC.

As shown in Fig 12a, there is good agreement for all sets until an altitude of 300 m. However, the greatest wind speeds were found in sets 3 and 4, and the best matches to the SODAR data (about 6 m s^-1) were found in sets 5 and 6. In Fig. 12b, regarding the passage of the sea breeze, the SODAR data showed a decrease in wind speed at an altitude of 200 m, with a minimum value (about 4.5 m s^-1) of approximately 75 m. Sets 5 and 6 also accompanied this slowing trend, with set 5 featuring more approximate values. The other sets showed no major changes compared with the previous schedule. Figure 12c shows that the wind speed from the SODAR data remained approximately constant (4.5 m s^-1) until an altitude of 200 m, then increased slightly (about 6 m s^-1) to approximately 400 m. Only sets 3 and 4 did not follow this trend of decreasing wind speed with increasing altitude.

Figure 13a shows that all sets matched the SODAR data between 200° and 215° in the wind direction until approximately 400 m. Figure 13b shows that the SODAR data for wind speed had a sudden change in wind direction between 75 and 125 m, taking the N-W direction below 75 m. Between 125 and 400 m, the direction was maintained as S-W, and above 400 m, it was in the S-E direction. Only set 6 followed this trend, especially at low altitudes, managing to capture the changes in a more satisfactory way. Figure 13c shows that SODAR data trended toward the N-W direction, but at higher altitudes. The sets that come closest are 2, 6, 1, and 5. Sets 3 and 4 did not change from the S-W direction.

These results reinforce the idea that sets 3 and 4 did not perceive the sea-breeze input, presenting a time delay compared with the frontal passage measured by the instruments. Until now, it was possible to verify that set 6 detected the sea-breeze input more efficiently, particularly according to the SODAR data and measurement of TIBL height with LIDAR data. Therefore, the next section shows near-surface horizontal wind field simulations considering only set 6.

3.4Horizontal wind field simulations

Figure 14 shows near-surface wind vectors from grids D3 (297 km × 297 km) and D4 (63 km × 81 km) based on the average wind directions estimated by the simulation by using set 6 at 1100, 1300, and 1500 UTC. The line indicates the sea-land interface.

Figure 14a shows a southerly surface wind flow at 1100 UTC and a high pressure system on the west side, showing less intensity next to the Dunkerque region. Figure 14b also shows a southerly surface wind flow at 1100 UTC and the pressure system next to the Dunkerque region. Figure 14c shows a high pressure system in the region close to Dunkerque, which allowed the acceleration of the winds into the region. Figure 14d shows a more intense influence of the high pressure system with an effective sea-breeze input. Figure 14e also shows the high pressure system over the sea in the region near Dunkerque, to maintain the sea-breeze input. However, this occurred with the opposition of the synoptic wind coming from the continent, as shown in Fig. 14f.

3.5Statistical analysis

Typical statistical measures (Hanna, 1989; Seigneur et al., 2000; Han et al., 2008), such as mean bias (MB), normalised mean bias (NMB), root mean square error (RMSE), normalised mean square error (NMSE), and correlation coefficient (R), were calculated to evaluate the accuracy of the predictions made by the WRF model. The best value for MB, NMB, RMSE, and NMSE is 0, and the best value for R is 1 (a negative R value means that as one variable increases in value, the second variable decreases in value).

Table 2 contains the performance statistics for the WRF 15-min averaged temperature at 2 m, wind speed at 10 m, and wind direction at 10 m, compared with the MSS data. The confidence limits for the performance measures were calculated by using Hanna (1989) boot-strap resampling software, which shows that the differences in the performance measures between the models are significant at the 95% confidence level.

The MB has an absolute value of less than approximately 0.4 m s^-1 for wind speed, 35° for wind direction, and 3.2℃ for temperature. The RMSE is approximately 0.7-1.0 m s^-1 for wind speed, 107-135° for wind direction, and 3.0-5.0℃ for temperature. Sets 3 and 4 tended to slightly overestimate the wind speed (positive mean bias), but WRF sets 1 and 2 had only a small mean bias. Sets 4, 3, and 2 had slightly lower RMSEs than sets 1, 5, and 6 regarding wind direction simulations (a lower value of 107° for set 4 vs a higher value of 135° for set 6). However, in such light wind situations, wind directions are known to be variable and unreliable. Note that these means were calculated by using scalar wind directions rather than vector wind directions. The sensitivity simulation for set 2 yields better statistics for temperature, even better than those of the other scenarios, except for the MB and NMB. For temperature, the R correlations ranged from 0.77 (set 4) to 0.89 (sets 2 and 6). For wind speed, the R correlations ranged from 0.66 (set 3) to 0.82 (set 5). For wind direction, the R correlations ranged from -0.35 (set 6) to 0.17 (set 2). The lower R for the wind direction by set 6 (-0.35) is due to the poor performance after 1800 UTC.

The poorer skill for the wind fields is mainly due to the inability of the currently used grid spacing to capture subgrid scale fluctuations induced by local topography and the underlying surface. Smaller correlations for the wind direction could also be associated with the limitation of the comparison method. Wind direction is actually a vector, whereas the comparison is intended for scalar values; thus, the correlations between the directions might not be reasonably reflected, especially when the wind direction is approximately 0° or 360°. The computation of statistical parameters is straightforward for wind speed, but the circular nature of wind direction makes it difficult to obtain the corresponding statistics. To avoid this we also used a modified wind direction, whereby 360° was either added to or subtracted from the predicted value to minimize the absolute difference between the observed and predicted wind directions (for details see the work of Lee and Fernando, 2004).

4Discussion

A good review about numerical studies on the sea breezes can be found in the work of Crosman and Horel (2010). It states that all studies reveal that higher values of land-surface sensible heat flux are associated with deeper and stronger sea breezes. At this point it is important to mention that the magnitude of the land-surface sensible heat flux influences all aspects of the circulation, constituting the fundamental driver of sea breezes (Crosman and Horel, 2010). Increasing the land-surface sensible heat flux tends to increase the intensity of sea breezes, which acts to increase the inland penetration distance. Nevertheless, as the landsurface sensible heat flux increases, turbulent convection also increases, which acts to annul the thermal gradient along the sea-breeze front. This process is known as turbulent frontolysis and decreases the inland penetration of the sea breeze front through a weakening of the horizontal temperature gradient (Ogawa et al., 2003).

In the past, a large number of simulations of sea breezes were realized with horizontal grid spacing greater than 2 km and lower-order ABL turbulence parameterizations. The largest deficiency of many early numerical models of sea breezes may be related to the unsuitable use of ABL turbulence (Yang, 1991), and may be the reason that most simulations are unable to represent the observed afternoon slowing of the sea-breeze front associated with turbulent frontolysis (Simpson et al., 1977). However, increasing the horizontal resolution in the simulations or using more sophisticated ABL parameterization schemes does not ensure improved results (Srinivas et al., 2007). Generally, the structures of sea breezes are well-represented, but the minor scale features and interactions between the sea breeze and the others variables are normally not captured.

The interaction among large-scale synoptic winds and sea breeze has been well reported in the literature e.g., (Zhong and Takle, 1993) by numerical-theoretical (Estoque, 1962; Pielke, 1974; Arritt, 1993; among others) and observational (Gilliam et al., 2004; among others) studies. Zhong and Takle (1993) showed that an onshore large-scale flow weakens the sea breeze. The impact of large-scale flows is greater than the thermal gradient caused by heating during the daytime. Furthermore, the wind direction of large-scale flows plays an important role in the evolution of sea breezes. Previous studies have shown that a perpendicular offshore large-scale wind strengthens sea-breeze perturbation by compressing the horizontal land-sea temperature gradient. Onshore geostrophic winds disturb this thermal gradient and consequently weaken the sea breezes (Azorin-Molina and Chen, 2008).

The ABL height (h) is important in atmospheric models because this height is used for other physical parameterizations. Explicitly, the methods for determining the ABL height are not identical to the designs tested. In the YSU design for unstable conditions, the height is determined as the first neutral level by checking the bulk Richardson number between the lowest level (z₁) and upper levels (Hong et al., 2006). In the ACM2 design, the method is similar to that of YSU; ABL height is detected as the height of the neutral floatation level where the bulk Richardson number for entrainment exceeds a critical value (Pleim, 2007b). Diffusive profiles are limited to heights below z=h in both designs because of the profiles predicted for the diffusion coefficient K. Thus, ABL height and temperature profiles are directly connected to these two designs. However, in the MYJ design, ABL height is detected as the height where the TKE reaches a sufficiently small value. Then, in this design, there was no direct connection between ABL height and temperature profile.

5Conclusions

The measurements obtained during the experiment and the model results identify changes in wind speed and direction, temperature decline, increases in humidity and an increase in vertical ascent in the vicinity of the sea-breeze front. Except for set 6, all simulations showed a delay in the sea breeze compared with the SODAR and LIDAR measurements. Both instrument recordings and model simulations identified the presence of the TIBL, which decreased continuously until its disappearance in the evening. The sea breeze had, during its entry, a synoptic wind opposition that attenuated the intensity and penetration. This provided speed and wind direction changes in a short period of time, with particular attention to soil temperature and neighboring regions that could not be detected by WRF simulations with the computational grid used (1-km grid cell resolution).

The local-closure MYJ scheme in combination with the land surface Noah scheme produced the best estimation of the ABL height in the presence of sea breeze in this study, and it captured the TIBL phenomena. The non-local closure ACM2 scheme presented the greatest delay compared with the YSU (non-local closure) and MYJ schemes in the estimation of sea breeze. In this work, the use of complex turbulence parameterizations (ACM2) did not improve the simulated mean and turbulent properties in the lower atmosphere. The results suggest that variations among the designs for determining breeze entry are mainly because of differences between the intensities of vertical mixing due to the ability of the soil parameterization to account for surface temperature. The surface processes, such as radiative exchanges and evapotranspiration, govern the partitioning of net radiation into sensible, latent and soil heat fluxes, which in turn strongly influence ground level air temperature and humidity, as well as ABL development. In particular, the prediction of the soil surface temperature and moisture content is critical for obtaining successful forecasts of heat and moisture exchange between the surface and the atmosphere. These results show that there is a need for thorough verification of the schemes using additional experimental observations, which are planned for the region.