Effect of Land-Use Change on the Urban Heat Island in the Fukuoka – Kitakyushu Metropolitan Area , Japan

In coastal cities, the effect of the sea breeze in mitigating the urban heat island (UHI) phenomenon has attracted attention. This study targeted the Fukuoka–Kitakyushu metropolitan area, the fourth largest metropolitan area in Japan which is also coastal. Doppler Light Detection And Ranging (LiDAR) observations were conducted in the summer of 2015 to clarify the transition of the wind field over the targeted area. To investigate the effects on the UHI of land-use change related to urbanization, the National Land Numerical Information (NLNI) land-use datasets for Japan in 1976 (NLNI-76) and 2009 (NLNI-09) were used in the Weather Research and Forecasting (WRF) model. The results of the simulation showed that most of the northern part of the Kyushu region became warmer, with an average increase of +0.236 ◦C for the whole simulation period. Comparing the two simulations and the Doppler LiDAR observations, the simulation results with the NLNI-09 dataset (for the year closest to the study period in 2015) showed closer conformity with the observations. The results of the simulation using NLNI-76 showed faster sea breeze penetration and higher wind velocity than the observations. These results suggest that the land-use change related to urbanization weakened the sea breeze penetration in this area.


Introduction
Japan is a mountainous island nation most of whose large cities are located in coastal areas and many of which are facing problems caused by urbanization.The urban heat island (UHI) phenomenon is a typical environmental problem encountered in dense urban areas in summer.The effect of the sea breeze in mitigating the UHI phenomenon has attracted attention in coastal cities, making further investigation of the environment of coastal urban areas important.
This study targeted the Fukuoka-Kitakyushu metropolitan area for investigation.This is the fourth largest metropolitan area after Tokyo, Osaka, and Nagoya, all of which are also coastal cities.Within these areas, sea breezes are important factors that mitigate the air-temperature rise during the summer.However, ongoing urbanization could be changing not only the mechanism of the energy balance on urban surfaces, but also the sea breeze system in large coastal cities.
The Fukuoka-Kitakyushu metropolitan area comprises the cities of Fukuoka and Kitakyushu in the Fukuoka prefecture located in the Kyushu region.The heat island monitoring report [1] indicated that the rise in annual average air temperature statistically derived with the linear least squares method from 1931 to 2014 in Fukuoka, in the central part of the targeted area, at 3.1 • C per century, was the second largest after that in Tokyo, at 3.2 • C per century.Moreover, the rise in minimum summer air temperature statistically derived with the linear least squares method from 1931 to 2014 in Fukuoka was the largest, at 3.9 • C per century.
To clarify the effects of land-use change related to urbanization in the Fukuoka-Kitakyushu metropolitan area, a meso-scale meteorological model was utilized.The term "meso" means intermediate between "macro" and "micro", and in the meteorological field, meso-scale means from a few kilometers to several hundred kilometers.To analyze the UHI in terms of meso-scale phenomena, e.g., motion of the atmosphere, heat transfer between land surface and atmosphere, solar radiation and infrared radiation and so on, meso-scale meteorological models are utilized.The Weather Research and Forecasting model (WRF) [2] developed by the National Center for Atmospheric Research (NCAR) in Boulder, CO, USA, is a meso-scale meteorological model that is widely utilized to analyze UHI phenomena.For example, Holt and Pullen analyzed the UHI in the New York City metropolitan area, USA, using the WRF and the Coupled Ocean-Atmosphere Mesoscale Prediction System (COMPAS) [3].Miao et al. analyzed the UHI and boundary layer structures in Beijing, China, with the WRF [4].Giannaros et al. investigated the spatial and temporal structure of the UHI in Athens, Greece, with the WRF [5].In Japan, Kusaka et al. utilized the WRF to evaluate the effect of future global warming in the Tokyo, Osaka, and Nagoya metropolitan areas [6].Hisada et al. investigated the relationship between the characteristics of sea breeze and land use in the Fukuoka metropolitan area using the WRF [7].To evaluate the effect of land-use change, the WRF is also adopted in this study.
The National Land Numerical Information (hereafter referred to as NLNI) land-use datasets provided by the National Land Information Division, National Spatial Planning and Regional Policy Bureau, Ministry of Land, Infrastructure, Transport and Tourism (MLIT) of Japan [8] were used to represent the land-use change in the Fukuoka-Kitakyushu metropolitan area, using a Land Surface Model (LSM) with the mosaic option in the WRF.The NLNI land-use datasets were published for several years, with the earliest dataset published in 1976 and the latest published in 2009.These were adopted to represent the land-use change over about three decades.The investigation period was set in the summer of 2015.To clarify the effect of land-use change, the same meteorological and simulation conditions were applied to the two simulation cases, except for the land-use fraction dataset.
Additionally, upper air observations using Doppler Light Detection And Ranging (LiDAR) were conducted to clarify the transition of the wind field over the Fukuoka-Kitakyushu metropolitan area.The Doppler LiDAR system emits laser pulses to the atmosphere and detects the reflected light back-scattered from aerosols.Wind velocity along the line of sight is then obtained using the Doppler frequency shifts from the back-scattered light.In the field of urban climate, Doppler LiDAR systems are utilized to observe wind profiles over urban areas.For example, Lane et al. [9] and Drew et al. [10] observed wind profiles over central London, UK.In the Tokyo metropolitan area, some research on sea breezes flowing over the urban area has been conducted using upper air observations.For example, Yoshikado and Kondo carried out large-scale upper air observations using radiosondes and pilot balloons, to clarify the transportation of atmospheric pollutants affected by the sea breeze from Tokyo Bay in the 1980s [11].Following that study, the present author conducted upper air observations using two sets of Doppler LiDAR in 2011 [12].One instrument was set at approximately 12 km from the coastline of Tokyo Bay, and the other was set at approximately 34 km from the coastline.These sets of Doppler LiDAR observations captured the sea breeze penetration from Tokyo Bay.In 2012, Yoshikado et al. [13] carried out upper air observations using Global Positioning System (GPS) sondes and clarified the thermal structure over the Tokyo metropolitan area in summer.In the Fukuoka-Kitakyushu metropolitan area, this kind of upper air observation that focuses on the urban climate has been limited.For example, Hisada et al. [14] observed wind fields over Fukuoka city using Doppler SODAR (SOnic Detection And Ranging); however, compared with the Tokyo metropolitan area, the accumulation of observation data is not sufficient.Therefore, to clarify the transition of the wind field over the Fukuoka-Kitakyushu metropolitan area, Doppler LiDAR observations were conducted simultaneously with the numerical simulations.

Outline of the Meso-Scale Meteorological Model
The WRF ver.3.7.1 was used to analyze the atmospheric conditions of the targeted area, and a mosaic land-use categories option for the Noah LSM was incorporated from the WRF ver.3.6 [15].In the previous LSM, the land surface property of each analysis cell was defined by the dominant land-use category.If one analysis cell consisted of 60% "urban" land-use and 40% "forest" land-use, that cell was defined as "urban", and the characteristics of a forest were neglected.However, the mosaic LSM enables the representation of different land-use fractions in each analysis cell.In an already urbanized area, the amount of further conversion from other land uses to "urban" within a short time period may be small.Therefore, the mosaic LSM is useful in analyzing urban climate transitions caused by land-use change over several decades, which is not represented by the dominant land-use.

Outline of the National Land Numerical Information Land-Use Datasets
In this study, the NLNI land-use datasets were used to represent the progress of land-use change.
To make the NLNI land-use datasets suitable for the WRF LSM, the WRF Preprocessing System (WPS) was modified [16].The NLNI land-use datasets were published in 1976, 1987, 1991, 1997, 2006 and 2009.From these, the datasets that were published in 1976 (earliest) and 2009 (latest) were used in this study.The NLNI land-use datasets are stored on the basis of the Basic Grid Square (Third Area Partition) with a resolution of approximately 1 km 2 and a grid spacing with latitude and longitude of 30" and 45", respectively.The Basic Grid Square is widely used for statistics regarding population, land-use, and urban planning, among others, by national and local governments in Japan.A resolution of approximately 1 km 2 is suitable as input data for the meso-scale model.Each Basic Grid Square has fractional land-use data with different categories by year.Since each dataset comprised different land-use categories, the original ones were re-grouped into six categories of land-use to adapt to the WRF, as shown in Table 1.In the WRF model, the land-use dataset provided by the U.S. Geological Survey (USGS) (Reston, VA, USA) is used as the default.Figure 1a shows the geography of the Kyushu region, the main target of this study.Figure 1b shows the fraction of land in the "urban" category (shaded in Table 1) in the Fukuoka prefecture in 1976.The resolution of the figure is 1 km, and the degree of shading indicates the fraction of "urban" land-use within each 1 km 2 cell.Figure 1c shows the increment in the "urban" land-use fraction over time."Urban" land-use fractions in 1976 were subtracted from those in 2009; so that a positive value (red in Figure 1c) indicates an increase in "urban" land-use in that cell.Overall, the results obtained over the three decades did not show extremely rapid land-use change in the Fukuoka-Kitakyushu metropolitan area.In particular, the increments in "urban" land-use in the central part of Fukuoka city were very small, because this area had already been extensively urbanized by the 1970s.However, the "urban" land-use area had increased in the surrounding areas due to urban expansion.land-use within each 1 km 2 cell.Figure 1c shows the increment in the "urban" land-use fraction over time."Urban" land-use fractions in 1976 were subtracted from those in 2009; so that a positive value (red in Figure 1c) indicates an increase in "urban" land-use in that cell.Overall, the results obtained over the three decades did not show extremely rapid land-use change in the Fukuoka-Kitakyushu metropolitan area.In particular, the increments in "urban" land-use in the central part of Fukuoka city were very small, because this area had already been extensively urbanized by the 1970s.However, the "urban" land-use area had increased in the surrounding areas due to urban expansion.

Analysis Conditions
Figure 2 illustrates the analysis domains, which covered an area of 3000 km (east-west) × 2500 km (south-north) with a horizontal resolution of 25 km.The entire analysis domain was two-way nested into two-fold sub-domains: domain 2 (750 × 750 km with a resolution of 5 km) and domain 3 (200 × 180 km with a resolution of 1 km).The vertical domain was divided into 48 unequally spaced grids.The top was set by pressure at 5000 Pa (approximately 20 km in height), and the lowest grid height was approximately 25 m.
The summer of 2015 was selected as the study period.The simulations were started on 24 July, and the time integration was performed for 69 days.The first 15 days were for spin-up, and the results of the hourly output data for 54 days, from 8 August to 30 September 2015, were used in this study.This evaluation period covered the upper air observation period.The initial conditions were obtained from the National Centers for Environmental Prediction Final operational global analysis data (NCEP FNL), with a resolution of 0.25°.Relaxation lateral boundary conditions were applied to domain 1, and nest lateral boundary conditions were applied to domains 2 and 3.The surface temperature and the sea surface temperature, for setting the lower boundary conditions, were obtained from the NCEP FNL data.In this study, the following schemes were applied for all domains and all cases: microphysics: Thompson scheme [17]; long-wave radiation: Rapid Radiative Transfer Model (RRTM) scheme [18]; short-wave radiation: Goddard scheme [19]; surface layer: Eta similarity scheme based on Monin-Obukhov similarity [20]; land surface: Noah LSM model [21] with mosaic option [16]; and planetary boundary layer: Mellor-Yamada-Janjic scheme [20].Grell 3D ensemble cumulus parameterization [22] was applied only for domain 1.Urban canopy models were not applied for all

Analysis Conditions
Figure 2 illustrates the analysis domains, which covered an area of 3000 km (east-west) × 2500 km (south-north) with a horizontal resolution of 25 km.The entire analysis domain was two-way nested into two-fold sub-domains: domain 2 (750 × 750 km with a resolution of 5 km) and domain 3 (200 × 180 km with a resolution of 1 km).The vertical domain was divided into 48 unequally spaced grids.The top was set by pressure at 5000 Pa (approximately 20 km in height), and the lowest grid height was approximately 25 m.
The summer of 2015 was selected as the study period.The simulations were started on 24 July, and the time integration was performed for 69 days.The first 15 days were for spin-up, and the results of the hourly output data for 54 days, from 8 August to 30 September 2015, were used in this study.This evaluation period covered the upper air observation period.The initial conditions were obtained from the National Centers for Environmental Prediction Final operational global analysis data (NCEP FNL), with a resolution of 0.25 • .Relaxation lateral boundary conditions were applied to domain 1, and nest lateral boundary conditions were applied to domains 2 and 3.The surface temperature and the sea surface temperature, for setting the lower boundary conditions, were obtained from the NCEP FNL data.In this study, the following schemes were applied for all domains and all cases: microphysics: Thompson scheme [17]; long-wave radiation: Rapid Radiative Transfer Model (RRTM) scheme [18]; short-wave radiation: Goddard scheme [19]; surface layer: Eta similarity scheme based on Monin-Obukhov similarity [20]; land surface: Noah LSM model [21] with mosaic option [16]; and planetary boundary layer: Mellor-Yamada-Janjic scheme [20].Grell 3D ensemble cumulus parameterization [22] was applied only for domain 1.Urban canopy models were not applied for all domains because of the lack of detailed building parameters.Therefore, the effect of urbanization was represented by urban land-use, and the transition of urban morphology could not be considered.
domains because of the lack of detailed building parameters.Therefore, the effect of urbanization was represented by urban land-use, and the transition of urban morphology could not be considered.To evaluate the climatological effect of land-use changes over the three decades, numerical simulations were performed using two different land-use datasets.For simplicity, these two datasets are abbreviated as NLNI-76 for the National Land Numerical Information land-use dataset for 1976, and NLNI-09 for 2009.In order to evaluate the influence of the land-use datasets, the same meteorological data were adopted as the initial conditions for both cases.Furthermore, to evaluate the effects of the mosaic LSM and NLNI land-use dataset, a USGS land-use dataset, used in the WRF by default, was also adopted.For the USGS dataset only, the Noah LSM model without the mosaic option was adopted, with the other settings remaining the same as for NLNI-76 and NLNI-09.

Outline of the Doppler LiDAR System
For the observations, the Leosphere WINDCUBE V2 Doppler LiDAR system was set up at the Ohashi campus of Kyushu University (33.5597°N 130.4308° E), approximately 6 km from the coastline, with a height of 8 m above ground level (AGL).The location of the Doppler LiDAR observation site is shown as no.15 in Figure 1b,c.
The Leosphere WINDCUBE V2 Doppler LiDAR system emits eye-safe pulsed laser beams at fixed angles; four beams with zenith angles of 28° facing to the north, east, south and west respectively, and one vertical beam.The wind speed range is from 0 to 60 m/s and the resolution is 0.1 m/s.The wind direction resolution is 2°.The vertical resolution of the Doppler LiDAR is 20 m, and the range is from 40 m to 260 m.The data sampling rate is 1 Hz.In this study, measurements were averaged over 10-minute intervals.

Comparing Observation and Simulation Results of the Surface Air Temperature
The NLNI-09 utilized the latest land-use dataset, from the year (2009) closest to the simulation period in 2015.Therefore, to validate the simulation performance, the results of the NLNI-09 simulation were compared first with the observation data provided by the Japan Meteorological To evaluate the climatological effect of land-use changes over the three decades, numerical simulations were performed using two different land-use datasets.For simplicity, these two datasets are abbreviated as NLNI-76 for the National Land Numerical Information land-use dataset for 1976, and NLNI-09 for 2009.In order to evaluate the influence of the land-use datasets, the same meteorological data were adopted as the initial conditions for both cases.Furthermore, to evaluate the effects of the mosaic LSM and NLNI land-use dataset, a USGS land-use dataset, used in the WRF by default, was also adopted.For the USGS dataset only, the Noah LSM model without the mosaic option was adopted, with the other settings remaining the same as for NLNI-76 and NLNI-09.

Outline of the Doppler LiDAR System
For the observations, the Leosphere WINDCUBE V2 Doppler LiDAR system was set up at the Ohashi campus of Kyushu University (33.5597 • N 130.4308 • E), approximately 6 km from the coastline, with a height of 8 m above ground level (AGL).The location of the Doppler LiDAR observation site is shown as no.15 in Figure 1b,c.
The Leosphere WINDCUBE V2 Doppler LiDAR system emits eye-safe pulsed laser beams at fixed angles; four beams with zenith angles of 28 • facing to the north, east, south and west respectively, and one vertical beam.The wind speed range is from 0 to 60 m/s and the resolution is 0.1 m/s.The wind direction resolution is 2 • .The vertical resolution of the Doppler LiDAR is 20 m, and the range is from 40 m to 260 m.The data sampling rate is 1 Hz.In this study, measurements were averaged over 10-min intervals.

Comparing Observation and Simulation Results of the Surface Air Temperature
The NLNI-09 utilized the latest land-use dataset, from the year (2009) closest to the simulation period in 2015.Therefore, to validate the simulation performance, the results of the NLNI-09 simulation were compared first with the observation data provided by the Japan Meteorological Agency (JMA).To evaluate the effects of the mosaic LSM and the different land-use datasets, the results of the USGS case were also compared.In Fukuoka prefecture, the surface air temperatures at 14 sites were observed.The locations of the observation sites are shown in Figure 1b,c.Table 2 shows the mean error (ME) and root mean square error (RMSE) for the observed surface air temperatures and simulated air temperatures at 2 m AGL.The data for evaluation were hourly data for 54 days from 8 August to 30 September 2015.The domain used for this evaluation was domain 3, with a resolution of 1 km.The ME of the NLNI-09 simulation for the 14 sites was slightly positive (+0.128), while the ME of the USGS simulation for the 14 sites was negative (−0.731).These differences suggest that the USGS land-use dataset was unable to represent the urban area in the targeted area correctly.Table 3 shows the land-use fractions in the NLNI-09 dataset, and the dominant land-use categories in the USGS dataset, for the 14 grid cells with a resolution of 1 km that include the JMA observation sites.Comparing the NLNI-09 and USGS simulations, the dominant land-use categories in the NLNI-09 and the USGS datasets showed non-conformity.The USGS simulation underestimated the area of urban land-use.
The RMSE of the NLNI-09 simulation ranged from 2.553 to 3.473, and the average of the 14 sites was 3.247.The RMSE of the USGS simulation ranged from 2.979 to 3.748, and the average of the 14 sites was 3.391.The NLNI-09 simulation showed better results; however, the difference was small.These results suggest that significant errors occurred in the numerical simulation.Further investigation of these simulation errors is needed.
Hereafter, the results of the NLNI-76 and NLNI-09 simulations are discussed, to consider the effect of land-use change on the urban climate of the Fukuoka-Kitakyushu metropolitan area.

Surface Air Temperature Change Caused by Land-Use Change
Figure 3 shows the difference in air temperature at 2 m AGL between the NLNI-76 and NLNI-09 simulations.The hourly air temperature data for the 54 days from 8 August to 30 September 2015 were averaged for the land surface only.The average air temperatures for the NLNI-76 simulation were subtracted from those of the NLNI-09 simulation.Negative values (blue tones in Figure 3) indicated that the air temperature in NLNI-76 was higher than that in NLNI-09, and positive values (red tones in Figure 3) indicated that the air temperature in NLNI-76 was lower than that in NLNI-09.The RMSE of the NLNI-09 simulation ranged from 2.553 to 3.473, and the average of the 14 sites was 3.247.The RMSE of the USGS simulation ranged from 2.979 to 3.748, and the average of the 14 sites was 3.391.The NLNI-09 simulation showed better results; however, the difference was small.These results suggest that significant errors occurred in the numerical simulation.Further investigation of these simulation errors is needed.
Hereafter, the results of the NLNI-76 and NLNI-09 simulations are discussed, to consider the effect of land-use change on the urban climate of the Fukuoka-Kitakyushu metropolitan area.

Surface Air Temperature Change Caused by Land-Use Change
Figure 3 shows the difference in air temperature at 2 m AGL between the NLNI-76 and NLNI-09 simulations.The hourly air temperature data for the 54 days from 8 August to 30 September 2015 were averaged for the land surface only.The average air temperatures for the NLNI-76 simulation were subtracted from those of the NLNI-09 simulation.Negative values (blue tones in Figure 3) indicated that the air temperature in NLNI-76 was higher than that in NLNI-09, and positive values (red tones in Figure 3) indicated that the air temperature in NLNI-76 was lower than that in NLNI-09.
(a) Average for the whole period (from 01:00 to 24:00).As mentioned above, the central parts of Fukuoka and Kitakyushu (the darker colors in Figure 1b) were already urbanized in the 1970s, so the air temperature increases caused by land-use change As mentioned above, the central parts of Fukuoka and Kitakyushu (the darker colors in Figure 1b) were already urbanized in the 1970s, so the air temperature increases caused by land-use change in these areas have been small.However, because of urban expansion, the urban heat island phenomenon has expanded.Although some areas in Figure 3a show decreases in the air temperature, the air temperature difference averaged over the land surface of domain 3 for the whole period was an increase of +0.236 • C. Figure 3b,c shows the temperature rise in daytime and nighttime periods, in the same manner as Figure 3a.The temperature rise averaged over the land surface of domain 3 in the nighttime was +0.245 • C, which is larger than the daytime average of +0.225 • C.Moreover, the nighttime range, with a minimum value of −0.324 and a maximum value of +2.656 • C, is wide compared to the daytime range, which has a minimum value of −0.144 and a maximum value of +1.317 • C.

Transition of the Sea Breeze over the Urban Area
Sea breeze circulation is driven by the air pressure difference between the land and the sea surfaces, caused by the temperature difference between these surfaces.The temperature difference is caused by solar radiation.Therefore, the sea breeze circulation is clearly observed on clear days with a small synoptic scale pressure gradient.During the observation period, 11 September 2015 gave the best conditions to observe the sea breeze blowing across the urban area.On that day, the cloud cover was less than 1 okta for a whole day, and the pressure gradients were small.The gradient of mean sea level pressure between the Fukuoka district meteorological observatory and the Kumamoto local meteorological observatory was 0.7 hPa, with a distance of approximately 90 km in the north-south direction; and the gradient of mean sea level pressure between the Oita local meteorological observatory and the Nagasaki local meteorological observatory was 0.2 hPa, with a distance of approximately 170 km in the east-west direction.The locations of observatories are shown in Figure 1a.
The Doppler LiDAR observation site was located approximately 6 km from the coastline of Hakata Bay in a south-southeast direction.Therefore, in this study, the sea breeze was defined as a northerly wind ranging from the west-northwest to the north-northeast, with an average velocity of 2 m/s or greater based on the Doppler LiDAR observations at 12 heights from 40 m to 260 m with a vertical resolution of 20 m.
Figure 4 shows the time-height cross sections of the horizontal wind vectors of the Doppler LiDAR observations, and the numerical simulations with the NLNI-76 and the NLNI-09 land-use datasets, on 11 September 2015.The observation results were averaged over 10-minute intervals.Figure 5 shows vertical wind and potential temperature profiles at selected times.Here, the potential temperature is the temperature of an air parcel which changed adiabatically from its initial state to the standard pressure, 1000 hPa.Because the air temperature is affected by the air pressure, to investigate the thermal structures in the vertical direction, the potential temperature is an important parameter.in these areas have been small.However, because of urban expansion, the urban heat island phenomenon has expanded.Although some areas in Figure 3a show decreases in the air temperature, the air temperature difference averaged over the land surface of domain 3 for the whole period was an increase of +0.236 °C. Figure 3b,c shows the temperature rise in daytime and nighttime periods, in the same manner as Figure 3a.The temperature rise averaged over the land surface of domain 3 in the nighttime was +0.245 °C, which is larger than the daytime average of +0.225 °C.Moreover, the nighttime range, with a minimum value of −0.324 and a maximum value of +2.656 °C, is wide compared to the daytime range, which has a minimum value of −0.144 and a maximum value of +1.317 °C.

Transition of the Sea Breeze over the Urban Area
Sea breeze circulation is driven by the air pressure difference between the land and the sea surfaces, caused by the temperature difference between these surfaces.The temperature difference is caused by solar radiation.Therefore, the sea breeze circulation is clearly observed on clear days with a small synoptic scale pressure gradient.During the observation period, 11 September 2015 gave the best conditions to observe the sea breeze blowing across the urban area.On that day, the cloud cover was less than 1 okta for a whole day, and the pressure gradients were small.The gradient of mean sea level pressure between the Fukuoka district meteorological observatory and the Kumamoto local meteorological observatory was 0.7 hPa, with a distance of approximately 90 km in the north-south direction; and the gradient of mean sea level pressure between the Oita local meteorological observatory and the Nagasaki local meteorological observatory was 0.2 hPa, with a distance of approximately 170 km in the east-west direction.The locations of observatories are shown in Figure 1a.
The Doppler LiDAR observation site was located approximately 6 km from the coastline of Hakata Bay in a south-southeast direction.Therefore, in this study, the sea breeze was defined as a northerly wind ranging from the west-northwest to the north-northeast, with an average velocity of 2 m/s or greater based on the Doppler LiDAR observations at 12 heights from 40 m to 260 m with a vertical resolution of 20 m.
Figure 4 shows the time-height cross sections of the horizontal wind vectors of the Doppler LiDAR observations, and the numerical simulations with the NLNI-76 and the NLNI-09 land-use datasets, on 11 September 2015.The observation results were averaged over 10-minute intervals.Figure 5 shows vertical wind and potential temperature profiles at selected times.Here, the potential temperature is the temperature of an air parcel which changed adiabatically from its initial state to the standard pressure, 1000 hPa.Because the air temperature is affected by the air pressure, to investigate the thermal structures in the vertical direction, the potential temperature is an important parameter.Figure 4a,c shows that the Doppler LiDAR captured sea breeze penetration at 10:00 and the simulation results for the NLNI-09 dataset showed the same tendency, although the wind directions were not accurate.The results of the NLNI-76 predicted sea breeze penetration at 09:00 are given in Figure 4b. Figure 5a shows vertical profiles at 08:00, before the sea breeze blows into the Ohashi Doppler LiDAR observation site.Vertical wind profiles at 09:00 in Figure 5b clearly show the difference between the Doppler LiDAR, the NLNI-76 simulation, and the NLNI-09 simulation.At 10:00, in Figure 5c, the Doppler LiDAR captured the sea breeze penetration.At the same time, the vertical wind profile of the NLNI-09 simulation showed higher velocity.At 14:00, in Figure 5d, the simulation with NLNI-76 showed a potential temperature drop below 50 m AGL, and the simulation with NLNI-09 showed a higher wind velocity than NLNI-76.At 15:00, in Figure 5e, the sea breeze with the highest velocity was observed, and the NLNI-09 simulation showed good agreement with the observations.On the other hand, the NLNI-76 simulation showed a higher velocity than the Doppler LiDAR and NLNI-09 results.Moreover, the potential temperature profile of NLNI-76 showed an inversion below 150 m AGL.At 18:00, in Figure 5f, the potential temperature profiles showed a large difference between NLNI-76 and NLNI-09 below 100 m AGL.
Although the NLNI-09 simulation showed good agreement with the Doppler LiDAR observations at 15:00 (Figure 5e), generally, the numerical simulations overestimated the wind velocities near the land surface.Moreover, the wind directions of simulations were not accurate compared with the Doppler LiDAR observations.One possible reason is that an urban canopy model was not utilized in this study; therefore, the simulation results underestimated the friction caused by urban morphology.
Comparing the two simulation results, the NLNI-09 simulation, which is the more urbanized case, showed lower wind velocity and later sea breeze penetration than the NLNI-76 simulation.The results of this study suggest that the land-use change related to urbanization has delayed sea breeze penetration and weakened its wind velocity.However, this result partially contradicts previous researches.Yoshikado [23,24] investigated the relationship between UHI and sea breeze using a two-dimensional hydrostatic boundary-layer model, and summarized that the UHI delays the inland penetration of the sea breeze, and that in coastal regions the UHI intensifies the sea breeze velocity during its growing stage.Kusaka et al. [25] utilized a hydrostatic three-dimensional model to analyze land-use change in the Tokyo metropolitan area, and summarized that the progress of urbanization delays the sea breeze penetration from Tokyo Bay and makes the wind velocity of the sea breeze stronger in the inland area.Our previous study [26] in the Tokyo metropolitan area, using the meso-scale meteorological model MM5 [27], showed the delay of sea breeze penetration from Tokyo Bay to the inland area; however, the difference of wind velocity in the inland area was not significant.In this study, comparing NLNI-76 with NLNI-09, the increase in the UHI is shown to delay sea breeze penetration.This result is consistent with previous studies [23][24][25][26].Meanwhile, the increase in the UHI has weakened the velocity of the sea breeze.
In this study, the simulation results of the WRF overestimated the wind velocities near the land surface.The possible reasons of these errors are uncertainties about the analysis model and the initial conditions.In terms of the model uncertainty, as mentioned above, the urban canopy model was not utilized in this study, and it was possible for the analysis model to underestimate the friction caused by urban morphology.On the other hand, it was also possible for the uncertainty of the initial conditions to cause these errors.To utilize the meso-scale meteorological models such as the WRF, the initial conditions are obtained from global analysis data such as the NCEP FNL data; and these regular, gridded global analysis data are produced from discretely distributed observations on Earth.In this process, uncertainty in the initial conditions of the meso-scale meteorological models cannot be avoided.Lorenz showed that very small differences in the initial conditions will cause errors in weather prediction over time [28,29].One way of reducing the uncertainties in the initial conditions is an ensemble simulation; however, the ensemble simulation needs a number of simulations with perturbations of initial conditions.As with the initial conditions, a multi-model ensemble simulation is one way to reduce the uncertainties of the model.In this study, evaluating the effect of urbanization on sea breeze penetration quantitatively is difficult because of the uncertainties of the model and initial conditions; therefore, what are shown here are the results of qualitative comparison.Further investigations are needed, with the ensemble simulations using various conditions and various models for various urban areas, in order to clarify the effect of land-use change on sea breeze circulation.

Conclusions
This study investigated the Fukuoka-Kitakyushu metropolitan area, which is the fourth largest metropolitan area in Japan.Compared with the Tokyo metropolitan area, the capital of Japan and the largest metropolitan area in the world, the number of urban climate studies in the targeted area is limited.Therefore, to accumulate data on the urban climate in the Fukuoka-Kitakyushu metropolitan area, numerical simulations with a meso-scale meteorological model (the WRF) and upper air observations using Doppler LiDAR were carried out simultaneously in the summer of 2015.For the numerical simulation, the NLNI land-use fraction datasets for Japan were used to represent the progress of land-use changes, in place of the default USGS land-use dataset in the WRF.
For simulation results of the surface air temperature, the results using the NLNI land-use dataset for 2009 showed fewer errors, when compared to the observations by JMA, than those using the USGS land-use dataset.This comparison suggested that the USGS land-use dataset was unable to represent the urban area precisely.
To consider the effect of land-use changes over three decades on the urban climate of the targeted area, the NLNI land-use datasets for 1976 and 2009 were utilized in the WRF simulations with the mosaic LSM.Based on the surface air temperature, most of the northern part of the Kyushu region became warmer because of land-use change related to urbanization, with an average increase over the land surface of +0.236 • C for the whole simulation period.The temperature increases in the areas surrounding Fukuoka city and Kitakyushu city were significant because of the urban expansion.It is possible that the land-use changes in the Fukuoka-Kitakyushu metropolitan area also had an effect on the sea breeze penetration from Hakata Bay to Fukuoka city.Comparing the two simulation results and the Doppler LiDAR observations, the simulation results using the NLNI-09 dataset (for the year closest to the study period in 2015) showed greater agreement with the observations.The results of the NLNI-76 simulation showed faster sea breeze penetration and higher wind velocity compared with the observations.These results suggest that the land-use changes have weakened sea breeze penetration.However, this tendency partially contradicts other research.To expand the knowledge base for the mitigation of urban heat island phenomena by sea breezes in coastal cities in the summer, additional investigations are needed and it will be necessary to reduce the uncertainties of the simulations used.

( a )
Geography of Kyushu region and the locations of Japan Meteorological Agency (JMA) observatories.(b) "Urban" land-use fraction in 1976.

Figure 1 .
Figure 1.Geography of targeted area and "urban" land-use fraction in Fukuoka prefecture: (a) The geography of the Kyushu region.The locations of the Fukuoka district meteorological observatory, the Kumamoto local meteorological observatory, the Oita local meteorological observatory and the Nagasaki local meteorological observatory are also shown.(b) The "urban" land-use fraction in 1976.The degree of shading indicates the fraction of "urban" land-use within each grid cell.(c) Increase in the "urban" land-use fraction.The red indicates the degree of increase in the "urban" land-use fraction from 1976 to 2009, and the blue indicates the decrease.Numbers from 1 to 14 enclosed within circles indicate the observation sites of the Japan Meteorological Agency (Tokyo, Japan).Number 15 indicates the Doppler LiDAR observation site.

Figure 1 .
Figure 1.Geography of targeted area and "urban" land-use fraction in Fukuoka prefecture: (a) The geography of the Kyushu region.The locations of the Fukuoka district meteorological observatory, the Kumamoto local meteorological observatory, the Oita local meteorological observatory and the Nagasaki local meteorological observatory are also shown.(b) The "urban" land-use fraction in 1976.The degree of shading indicates the fraction of "urban" land-use within each grid cell.(c) Increase in the "urban" land-use fraction.The red indicates the degree of increase in the "urban" land-use fraction from 1976 to 2009, and the blue indicates the decrease.Numbers from 1 to 14 enclosed within circles indicate the observation sites of the Japan Meteorological Agency (Tokyo, Japan).Number 15 indicates the Doppler LiDAR observation site.

Figure 2 .
Figure 2. Analysis domains for the Weather Research and Forecasting model (WRF) simulation.Domain 1 covers all of Japan, domain 2 covers the entire Kyushu region, and domain 3 covers the northern part of the Kyushu region, which was the main target area in this study.

Figure 2 .
Figure 2. Analysis domains for the Weather Research and Forecasting model (WRF) simulation.Domain 1 covers all of Japan, domain 2 covers the entire Kyushu region, and domain 3 covers the northern part of the Kyushu region, which was the main target area in this study.

Figure 3 .
Figure 3. Differences in averaged air temperature between NLNI-76 and NLNI-09 over 54 days from 8 August to 30 September 2015.The blue tones show a temperature drop, and the red tones show a temperature rise from NLNI-76 to NLNI-09: (a) Difference in averaged air temperature for the whole period.The spatial average for this domain is +0.236 °C, with a minimum value of −0.171 °C and a maximum value of +2.049 °C.(b) Difference in averaged air temperature for the daytime period, from 09:00 to 18:00.The spatial average for this domain is +0.225 °C, with a minimum value of −0.144 °C and a maximum value of +1.317 °C.(c) Difference in averaged air temperature for the nighttime period, from 19:00 to 08:00 the next day.The spatial average for this domain is +0.245 °C, with a minimum value of −0.324 °C and a maximum value of +2.656 °C.

Figure 3 .
Figure 3. Differences in averaged air temperature between NLNI-76 and NLNI-09 over 54 days from 8 August to 30 September 2015.The blue tones show a temperature drop, and the red tones show a temperature rise from NLNI-76 to NLNI-09: (a) Difference in averaged air temperature for the whole period.The spatial average for this domain is +0.236 • C, with a minimum value of −0.171 • C and a maximum value of +2.049 • C. (b) Difference in averaged air temperature for the daytime period, from 09:00 to 18:00.The spatial average for this domain is +0.225 • C, with a minimum value of −0.144 • C and a maximum value of +1.317 • C. (c) Difference in averaged air temperature for the nighttime period, from 19:00 to 08:00 the next day.The spatial average for this domain is +0.245 • C, with a minimum value of −0.324 • C and a maximum value of +2.656 • C.

Figure 4 .
Figure 4. Time-height cross sections of horizontal wind vectors at the observation site on 11 September 2015: (a) Doppler LiDAR observation, (b) Numerical simulation with NLNI-76, and (c) Numerical simulation with NLNI-09.

Figure 4 .
Figure 4. Time-height cross sections of horizontal wind vectors at the observation site on 11 September 2015: (a) Doppler LiDAR observation, (b) Numerical simulation with NLNI-76, and (c) Numerical simulation with NLNI-09.

Figure 5 .
Figure 5. Vertical wind profiles from Doppler LiDAR observations and from simulations with NLNI-76 and NLNI-09, and potential temperature profiles from simulations with NLNI-76 and NLNI-09, on 11 September 2015: (a) At 08:00, before the sea breeze blows into the Ohashi Doppler LiDAR

Figure 5 .
Figure 5. Vertical wind profiles from Doppler LiDAR observations and from simulations with NLNI-76 and NLNI-09, and potential temperature profiles from simulations with NLNI-76 and NLNI-09, on 11 September 2015: (a) At 08:00, before the sea breeze blows into the Ohashi Doppler LiDAR observation site; (b) At 09:00, the simulation with NLNI-76 captured the sea breeze, while the Doppler LiDAR did not; (c) At 10:00, the Doppler LiDAR captured the sea breeze; (d) At 14:00, the simulation with NLNI-76 showed a potential temperature drop below 50 m above ground level (AGL); (e) At 15:00, the observed sea breeze became the strongest on that day; and (f) At 18:00, the potential temperature profile difference showed a large difference between NLNI-76 and NLNI-09 below 100 m AGL.

Table 1 .
Land-use categories in the Weather Research and Forecasting model (WRF) and National Land Numerical Information (NLNI) land-use datasets.

Table 2 .
Mean error (ME) and root mean square error (RMSE) for the surface air temperature when comparing the JMA observations and the simulation results for NLNI-09 and USGS.

Table 3 .
Land-use fractions in the NLNI-09 case, and the dominant land-use category in the USGS case, for the 14 grid cells with a resolution of 1 km that include the JMA observation sites.For the NLNI-09 case, the dominant land-use category is shaded.