Journal Archive

The Suita City government had selected two regions within the city where UHI mitigation strategies shall be implemented (Suita City, 2016). One of these regions is Esaka, the analysis area, which is a highly urbanized district with major roads, residential and industrial buildings. Prior to simulation, a 3D model of Esaka (4.758°N, 135.497°E) was created as shown in Fig. 1. The 3D model was then divided into meshes with a grid interval of 5 m. Fig. 2 shows the domain size and the analyzed area surrounded by the solid line. There were 22 buildings with the highest building being 56 m. The total mesh number in the CFD model was 80×77×47 including the virtual space. For the PT model, which is one dimensional, the total mesh number was 16 with a grid interval of 0.01 m.

Fig. 1. 
Esaka district (4.758°N, 135.497°E) as seen from Google Earth (2015) and its corresponding 3D model.

Fig. 2. 
Mesh view of the calculated domain.

The Weather Research Forecasting (WRF) model version 3.5.1 (Skamarock et al., 2008) was used to determine the boundary conditions of the CFD model. Fig. 3 shows two WRF domains: the coarser domain (D1) covering Kansai region with 90×90 grid cells and horizontal grid resolution of 3 km, and the finer domain (D2) covering Osaka Prefecture, which includes Suita City with 90×90 grid cells and horizontal grid resolution of 1 km. The vertical domain consists of 30 layers from the surface to 100 hPa with the middle height of the first, second and third layers are about 28 m, 92 m and 190 m above the surface, respectively. Except for the calculation period, the WRF configurations were essentially the same as those in Shimadera et al. (2015), the Yonsei University planetary boundary layer scheme (Hong et al., 2006), the WRF single-moment 6-class microphysics scheme (Hong and Lim, 2006), the Noah land surface model (Chen and Dudhia, 2001), the rapid radiative transfer model (Mlawer et al., 1997) for the long wave radiation, and the shortwave radiation scheme of Dudhia (1989) with initial and boundary conditions derived from the mesoscale model grid point value data by the Japan Meteorological Agency (JMA).

Fig. 3. 
The WRF domains showing Kansai region as D1 and Osaka Prefecture, which includes Suita City, as D2.

The WRF simulation was conducted using online one-way nesting in the two domains for August 8-9, 2011 with a 7-day initial spin-up. This period was chosen because there was clear and calm weather due to the high pressure system that formed above Japan (Fig. 4a). Because of the lack of meteorological observatory in Suita City, the WRF results in the study period were compared with observation data of JMA at the Osaka meteorological observatory at (135.518°E, 34.682°N) in Osaka City (Fig. 4b). WRF fairly well simulated temporal variations of air temperature and wind, including daytime stronger wind caused by prevailing southwesterly sea breeze and nocturnal calm condition. The WRF results at Esaka, Suita City for the CFD boundary conditions were generally similar to those at the Osaka meteorological observatory.

Fig. 4. 
Weather map at 0900 JST on August 8, 2011 obtained from JMA (a); time series of observed and WRF-simulated temperature and wind at Osaka meteorological observatory and WRF-simulated values at Esaka, Suita on August 8-9 (b).

Owing to the coarser vertical resolution of WRF model compared with CFD, Monin-Obukhov Similarity Theory (MOST) was applied. The MOST defines the vertical distribution of wind speed and air temperature using hourly WRF data at heights 28 m and 92 m with a roughness length of 0.1 m. Fig. 5 shows the lateral boundary conditions of CFD for air temperature and wind components u and v at each vertical level. Air temperature was warmer in the surface layer especially during noon. Wind speed was stronger during daytime and wind direction was generally south-west because of prevailing sea breeze.

Fig. 5. 
Lateral boundary conditions of CFD for air temperature and u, v wind components at each vertical level.

The CFD model was composed of the conservation equations for momentum, continuity and mass. The k-ε turbulence model was used for turbulence while Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm (Pantakar, 1980) was utilized for pressure corrections. The CFD model calculated for air temperature, wind speed, air humidity, downward longwave radiation and shortwave radiation and supplied these parameters to the PT model. The CFD calculation was conducted for 24 h from 0500 Japan Standard Time (JST) on August 8 to 0500 JST August 9, 2011 with a time step of 0.5 s. Linear interpolation was applied to the hourly data obtained from the WRF model in order to update the lateral boundaries of the CFD model every time step.

There were two surface energy budget model incorporated into the CFD model. One is the building envelope model or BEM (Eq. (1)), used to calculate heat transfer from building surface to atmosphere. For a detailed discussion of the BEM refer to Cortes et al. (2015).

Here ε is emissivity, σ is Stefan-Boltzmann constant (W/m² K⁴), T_bso is the temperature of outer building surface (°C), f_i,j is view factor of each surface element j from element i, T_j is temperature of surface j, α_bs is building surface albedo, S↓ is shortwave radiation (W/m²), ρ_a is air density (g/cm³), c_p is specific heat at constant pressure (J/kgK), c_{H_bs} is bulk heat transfer coefficient, u is wind speed (m/s), T_air is air temperature (°C) and Q is heat flux (W/m²). The method for calculating solar radiation and view factor can be obtained from Ikejima et al. (2011).

Second is the PT model which simulates heat transfer from ground surface to the atmosphere. The PT model was composed of conservation equations for heat and soil moisture. The parameters used were volumetric water content θ (m³/m³), evaporation efficiency β and matric potential ψ (m). The β was defined from the evaporation equation (Eq. (2)) for bare soil developed by Kondo et al. (1990). The van Genuchten-Mualem model shown in Eq. (3)-(4) was used to express the relationship between ψ and θ (van Genuchten, 1980; Mualem, 1976).

Here E is evaporation rate (kg/m²s), C_E is bulk coefficient of evaporation, q_s is specific humidity (g/kg) at ground surface temperature T_s (°C), q_a is specific humidity of air (g/kg), θ_sat is saturated volumetric water content, x and n are van Genuchten curve-fitting parameters, m=1- (1/n), K is the hydraulic conductivity (m/s) and K_sat is saturated hydraulic conductivity (m/s). In this study we were interested in determining the maximum performance of WRP hence the initial θ was set at saturation level (Table 1) and we did not consider other scenarios like varying initial θ. The ground surface energy balance is shown in Eq. (5) where the left terms represent radiation and right terms represent sensible heat flux, latent heat flux, ground heat flux respectively.

where α_g is ground albedo, L↓ is downward longwave radiation (W/m²), H is sensible heat flux (W/m²), l is latent heat of vaporization (J/kg) and G is heat flux into the soil (W/m²). The heat fluxes were solved using Eq. (6)-(8).

The bulk coefficient for the sensible heat flux, C_{H_gs}, was assumed to be equal to C_E based on Kondo et al. (1990). The initial values of parameters used in PT model are listed in Table 1. The PT model calculation time step was 60 s.

Fig. 6 shows the model integration and the variables calculated in each system. Using the coupled CFD-PT model, two cases of 24 h unsteady analysis were simulated. In case 1, asphalt was used as the pavement material of all ground surface. In case 2, WRP was used as the pavement material of the main street and asphalt for the rest of the ground surface.

Fig. 6. 
The calculation flow and coupling of meteorological models, CFD, BEM and PT model.

Examining the average diurnal variation in ground surface temperature within the analysis area (Fig. 7), it can be seen that T_s in case 2 was consistently lower throughout the day. As Cortes et al. (2016) pointed out, the cooling property of WRP compared with asphalt can be attributed to three factors namely water content, albedo and thermal conductivity. However, by creating a graphical representation of the diurnal variation in T_s, it was discovered that shadow effect also contributed to a decrease in surface temperature (Fig. 8). True for both case 1 and case 2, T_s on the western side of the buildings decreased at 0800 JST. During this time the shadow formed on the west as the sun rose in the east.

Fig. 7. 
The diurnal variation in ground surface temperature, T_s, within the analysis range for 24-h time period.

Fig. 8. 
The x-y view of ground surface temperature, T_s, in each case.

Here we focus the discussion to the difference between main street T_s in case 1 and case 2. Results show that T_s in case 2 main street decreased by 3.1°C at 0800 JST. We attribute this to the greater water content of WRP compared with asphalt which led to evaporative cooling and increased latent heat flux (discussed in section 3.3). This agrees with the study of Asaeda (1993) on evaporation in bare soil where the transport of water vapor inside soil affects subsurface distribution of temperature greatly. Due to evaporation, bare soil is cooler than covered surfaces such that increase in thickness of the covering material causes temperature increase and higher heat stored. Similar to bare soil, the ability of WRP to hold water and its porosity allow for water transport thus affecting subsurface distribution of temperature. It can be seen in Fig. 9 that at 0800 JST, the average water content of WRP was 0.08 which was still near to the saturation level of 0.09. At 1200 JST the water content decreased to 0.05, the greatest change throughout the day. This can be expected because at noontime, the sun was most intense. Increased solar radiation further promoted evaporation of water from surfaces. It also explains the greatest difference in T_s between case 1 and case 2; during this time, T_s in case 2 main street decreased by 13.8°C. From 1600 JST to 2000 JST the surface water content remained at 0.04 but the surface temperature continued to decrease. During this time, T_s of case 2 main street decreased by 5.7°C and 2.9°C. respectively. We attribute this decrease in T_s despite the constant surface water content to both higher albedo and lower thermal conductivity of WRP compared with asphalt (Table 1). Both factors allowed reflection of sunlight during daytime and minimized the release of heat during night time. These estimates of main street T_s between case 1 and case 2 were proximate to the observations of Cortes et al. (2016) on the behavior of asphalt and WRP in an outdoor set-up.

Fig. 9. 
Diurnal variation in water content (θ) in case 2 main street surface.

In this section we discuss the effect of WRP on diurnal variation in air temperature. The CFD-PT model estimated an overall decrease in air temperature at 1.5 m within the analysis area (Fig. 10). The average air temperature difference over main street was computed as case 1 T_air - case 2 T_air. Results show that air temperature in case 2 decreased by 0.05°C at 0800 JST, 0.3°C at 1200 JST, 0.1°C at 1600 JST and 0.06°C at 2000 JST. The air temperature decrease was most pronounced at 1200 JST which is proportional to the degree of surface temperature cooling. Looking at the x-y section in Fig. 11, the decrease is particularly noticeable on the western side of the main street area. We attribute this to the stronger wind towards north-west direction which is evident in the wind profile.

Fig. 10. 
The diurnal variation in air temperature at 1.5 m above ground surface.

Fig. 11. 
The wind profile and effect of WRP on air temperature.

By examining the x-z section, we can see greater air temperature decrease above WRP surface (main street) compared with asphalt surface (none main street areas) confirming the occurrence of evaporative cooling. Moreover, there was also vortex formation within the street canyon. This vortex allowed for mixing of air that eventually caused cooler air from the ground to rose.

Fig. 12 shows the diurnal variations in ground surface energy fluxes within the analysis area. Results show an overall increase in latent heat flux and decrease in sensible heat flux, conductive heat flux and upward longwave radiation throughout the day. Due to the significant decrease in main street T_s, a comparison of the surface energy fluxes in the same area was also carried out. The ground surface energy flux difference over main street was computed as case 1 heat flux- case 2 heat flux. The model estimated an average increase in latent heat flux by up to 51 W/m² at 0800 JST, 255 W/m² at 1200 JST, 68 W/m² at 1600 JST and 34 W/m² at 2000 JST. This diurnal variation in latent heat flux was also relative to the diurnal variation in water content discussed previously. As the sun started to become intense from morning to noon, evaporation rate also increased dramatically. Consequently, this increase in latent heat flux caused a particularly slight difference between T_s and T_air, allowing for a distinct decrease in sensible heat flux. As expected from the increase in latent heat flux, maximum decrease in sensible heat flux occurred at 1200 JST. During this time sensible heat flux in case 2 decreased by 465 W/m². The longwave radiation also decreased by 97 W/m².

Fig. 12. 
A comparison of ground surface energy fluxes between case 1 and case 2.