## Abstract

The atmospheric boundary exhibits an obvious diurnal cycle. The daily cycle of climate variation has a significant effect on urban airflow, and an understanding of it is very important for city-scale environmental control. A new and simple daily cycle temperature boundary condition for simulations of urban airflows with computational fluid dynamics (CFD) is described herein. An analytical surface temperature formula was obtained after simplifying longwave radiation and sensible heat flux terms. The formula provides a reasonably good prediction of ground surface temperatures on sunny days without adding much complexity. The accuracy of the prediction of daily surface temperature variations in homogeneous soils was evaluated with a benchmark experiment and with weather station data. The new boundary condition was implemented in the commercial CFD software Fluent for a city-scale model. The implementation demonstrates the importance of including diurnal temperature profiles and atmospheric boundary conditions in CFD simulations of urban plumes in which an ideal urban heat island circulation occurs.

## 1. Introduction

Understanding the detailed microscale environment in a multiscale environment (e.g., a city) has become a key topic of study in recent years (Martilli 2007). The general models used for city-scale wind flows can be divided into two groups. The first is based on the traditional mesoscale weather forecasting models, such as the Weather Research and Forecasting (WRF; Model Skamarock et al. 2005) or the fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5; Dudhia et al. 2000). By downscaling the mesoscale models or by coupling them indirectly with smaller-scale models, the mesoscale models can approximately capture the city-scale environment. The grid size of the downscaling mesoscale model may be insufficient for some applications such as studies of the wind flow in a street. For example, the high horizontal resolution of the downscaled mesoscale model simulated by WRF coupled with an urban canopy layer is around 500 m (Tewari et al. 2006). Nicholls (1993) used the Colorado State University Regional Atmospheric Modeling System to increase the resolution to 225 m, which is greatly dependent on the resolution of the land-use type and is not sufficiently accurate for specific buildings in a city.

The other group is based on computational fluid dynamics (CFD) models including Reynolds-averaged Navier–Stokes or large-eddy-simulation turbulence models, which are capable of solution at the fine grid scale (from 1 to100 m) (Xie 2011; Blocken and Gualtieri 2012; Martilli 2007; Wang 2009; Zhang et al. 2014). Although consideration of the 24-h daily temperature cycle is long established within the mesoscale models, most CFD simulations are carried out under neutral atmospheric conditions without considering the impact of atmospheric stability in cities (Lien and Yee 2004; Riddle et al. 2004). Some studies have considered stable or unstable atmospheric conditions, but none have focused on the evolution of atmospheric stability, which changes on a daily cycle. A full-cycle simulation of diurnal changes is also rare in the applications of CFD models because of the limited capability to simulate the city-scale environment. Kristóf et al. (2009) successfully proposed a method for upscaling a CFD model to simulate a city-scale environment. Wang and Li (2016) improved this model, which is referred to as a city-scale CFD (csCFD) model. The csCFD model is used here for further study for its effectiveness in handling the complex geometry that is a feature of urban heterogeneity and also for its ability to simulate the city-scale environment at a fine resolution (down to 1 m).

In reality, surface and air temperatures change every second, such as at the street level. The daily cycle of climate variation has a crucial impact on city-scale environmental control, especially in the planetary boundary layer (PBL). The PBL’s typical diurnal cycle is shown in Fig. 1. The importance of such diurnal temperature profiles has been revealed in numerous studies. The same pollutant source from a chimney disperses differently over the course of a day because of the differences in atmospheric stability and wind profiles (Arya 1999). The change in diurnal temperature profiles results from changes in the Earth surface energy balance—a topic that has been studied extensively and for which many models have been proposed (Piringer et al. 2002).

The energy balance process is very complicated because it involves all heat transfer modes, including radiation, convection, conduction, evapotranspiration, and so on. Each of these modes consists of several physical and chemical processes. Some terms in the surface energy balance budget are coupled with each other and with all of the atmosphere’s physical and chemical processes. The complete equation of the entire energy balance without any simplification is rarely used. The energy budget model has been simplified to varying degrees for different purposes. The aim of many surface energy balance models is to provide the surface temperature for mesoscale models such as WRF, MM5, and so on; therefore, the characteristics of the geometry of the buildings and the radiation exchange have been simplified and parameterized (Masson 2000; Kusaka and Kimura 2004; Oleson et al. 2008). With the advent of more powerful computers and improved understanding of the surface energy balance, the models were developed further, from the slab model (Wullschleger et al. 1991), to the single-layer model (Kusaka et al. 2001; Kanda et al. 2005), to the multilayer urban canopy models (Masson 2006). Vegetation (Lee and Park 2008; Kawai et al. 2009) and urban hydrological processes (Cheng and Wang 2002; Lemonsu et al. 2007) were recently added to the urban canopy model. From the microscale perspective, a few studies have focused on the microscale temporal distribution of the surface temperature in cities or building clusters. In these models, the three-dimensional geometry and radiation exchange are considered in a more detailed manner and are less simplified than those in mesoscale surface parameterizations (Krayenhoff and Voogt 2007; Yang and Li 2015).

Although it is necessary to include the variations in daily surface temperature, the existing CFD models or simulations do not consider the diurnal cycle of the atmospheric boundary layer (Blocken and Gualtieri 2012). Coupling the existing surface temperature model with CFD models is challenging. The commonly used thermal boundary conditions include constant temperature (Dirichlet boundary condition), constant heat flux (Neumann boundary condition), and the mixed type, which considers both convection and shortwave and longwave radiation from the sun and surrounding buildings. Implementation of an atmospheric boundary layer as boundary conditions is not entirely straightforward in a CFD simulation, mainly because of the incompatibility of the parameterizations used in mesoscale meteorological models and the boundary conditions required for a traditional CFD simulation (Toparlar et al. 2015). In both of the surface energy balance models used in the mesoscale models (e.g., Masson 2000) and in the dedicated three-dimensional surface temperature models (e.g., Yang and Li 2015), partial differential equations must be numerically solved. In comparison with simple analytical solutions, more computation time is required for iteration, which is not suitable when the surface temperature is changing with small time steps or when it is changing at many surfaces such as building walls and roofs. The computation time increases rapidly when the simulated time or surface increases. The purpose here is not to forecast the local- or city-scale environments as a general mesoscale model does but rather is to understand the parameters that govern daily climate variations. Therefore, what is needed is not a complex surface temperature model that governs through a system of partial differential equations that require many iterative steps to converge, but a simple and effective approach that will not substantially increase the computation time and will predict the surface temperature with sufficient accuracy. This can be accomplished by specifying the surface temperature and surface heat flux in CFD without the need for coupling with the ground heat transfer. This method is also useful for studies for which the computation time needs to be minimized. For example, public events need a fast response.

Williamson (2004) demonstrated the possibility of analytically solving the governing equation of the surface energy balance, but the longwave radiation term was simplified with a constant gradient of the linearized function. However, the sky temperature should be used instead of the air temperature in the longwave radiation calculation (Silva et al. 2009). In this study, the longwave radiation term was improved by including the sky temperature and improving the gradient term. Last, we derived a more accurate analytical solution, referred to as the simple daily cycle temperature (sDCT) boundary condition. It provides reasonably good results without adding much complexity. Note that the air temperature is only known at previous time steps in a time-dependent simulation; thus, there is a need to address the coupling of surface temperature with air temperature and the proper time step.

In this study, we chose to specify the sensible heat flux instead of the surface temperature as the boundary condition so that energy is easy to conserve in the calculations. We implemented the new boundary condition in the csCFD model (Wang and Li 2016). The csCFD model can treat the urban area as porous media and fully resolve some of the complex geometry (Wang et al. 2017), and the grid resolution is finer than the general mesoscale model. The feasibility of treating the urban area as porous media for building clusters with or without heat has already been demonstrated in the literature (Hang and Li 2010; Hang et al. 2012; Hu et al. 2012). Relative to a parameterized urban canopy model, using a porous model can provide a more precise wind environment and more precise air temperatures. As the first step, a boundary condition is developed to predict surface temperatures over a flat surface. It can be easily extended to slopes and to typical urban areas with buildings. The simple boundary condition can be used in combination with the three-dimensional (3D) radiation model of Yang and Li (2013) for the prediction of wind flows in building clusters or a city. We first aimed to evaluate the accuracy of the prediction of daily surface temperature variation in homogeneous soils. An example has been given to demonstrate the importance of including diurnal temperature profiles and diurnal atmospheric boundary conditions in csCFD simulations of urban plumes when the synoptic wind is weak.

## 2. Methods

We aim to improve the prediction of the urban surface temperature, which is an important boundary condition for the modeling of the atmospheric boundary layer. The surface discussed in this work is the same as the ground-level surface defined in Voogt and Oke (1997). The surface can be bare soil, concrete, grass, or other types, which differ in terms of their thermal properties. Assume that there is no heating source from Earth’s interior and that the only forcing in the interface between land and atmosphere is solar radiation. First, we considered sunny days when the influence of evaporation is insignificant. The effects of evaporation on engineered pavement can be ignored if there are at least two continuous sunny days (Wang et al. 2013). For a natural surface, sunny days are acceptable. The general surface energy balance can be expressed as

where *Q*_{NR} is the net radiation (W m^{−2}), *Q*_{SH} is the sensible heat flux (W m^{−2}), *Q*_{LE} is the latent heat flux (W m^{−2}), *Q*_{G} is the soil heat flux (W m^{−2}), and *Q*_{AN} is the anthropogenic heat (W m^{−2}). Both *Q*_{NR} and *Q*_{G} are defined as positive upward, and the other two are defined as positive downward. Hence, the net radiation *Q*_{NR} is balanced by the sum of the sensible heat flux into the air, the latent heat of evaporation, and the soil heat flux at the surface. Each term is described as shown in appendix A.

We rewrite the surface energy balance Eq. (1) as

The net radiation *Q*_{NR} is the sum of all of the shortwave and longwave radiation, which includes the incoming and outgoing shortwave radiation *Q*_{SW}, the incoming longwave radiation , and the outgoing longwave radiation :

The shortwave radiation can be calculated as

where *α*_{l} is the albedo of the surface (or land), which changes with different building morphologies, *α*_{a} is the albedo of the atmosphere, *C*_{A} is the fraction of the atmosphere, and *Q*_{s} is the solar radiation at the top of the atmosphere.

The solar radiation at the top of the atmosphere is related to location and time:

where *S*_{*} is the solar constant (=1368 W m^{−2}). In addition, *ω* = [(*t*/3600) − 12], *λ* is the local latitude, *δ* is the current declination of the sun, *t* is the local time (s), and = *π*/12. The shortwave radiation is an independent term.

By substituting Eqs. (3)–(5) and all of the simplifications or parameterizations in appendix A into Eq. (2), we obtain the following heat balance equation at the ground surface:

In this equation, the left-hand term is the effective heat storage of the ground surface and the right-hand terms are the solar heat gain, the net thermal radiation including sky radiation, the convective heat, and any possible anthropogenic heat generation, respectively. Latent heat from evaporation or condensation is ignored. The convective heat transfer coefficient *h*_{sh} is determined by the wind flow along the surface and by the surface temperature gradient, which can either be obtained with the CFD method or provided using an empirical formula; the latter approach is adopted here, as shown in appendix A.

Equation (6) can be rewritten as

where *A* = *h*_{l}*C*_{l}, *B* = *h*_{rad} + *h*_{sh}, *C* = (1 − *α*_{l})(1 − *α*_{a})(1 − *C*_{A}) cos*λ* cos*δS*_{*}, *D* = −*E* − *F*, *E* = (1 − *α*_{l})(1 − *α*_{a})(1 − *C*_{A}) sin*λ* sin*δS*_{*}, and *F* = (*h*_{rad}*γ* + *h*_{sh})*T*_{a} + *Q*_{AN}. The term *D* is initially assumed to be a constant to obtain the analytical solution, but in reality *D* is not a constant. The variables upon which the term *D* depends—including the air temperature *T*_{a}, the dewpoint temperature *T*_{dew}, and the reference velocity *U*_{10}, among others, as explained in appendix A—also change with time. Therefore, to further improve the accuracy of the analytical solution, the piecewise functions of *T*_{a}, *T*_{dew}, and *U*_{10,} and other variables are adopted here. The hourly results are used for the subsequent studies in this paper. If the results are obtained at a high temporal resolution, they could also be used here to improve accuracy.

For 24-h cycles, the shortwave radiation *Q*_{SW} can be expressed as

where *t*_{1} is the local time of sunrise and *t*_{2} is the time of sunset, calculated using

The current time is *t*. The values for *t*, *t*_{1}, and *t*_{2} are given in seconds. In period 1 (*t* ∈ [0, *t*_{1})) and period 3 (*t* ∈ (*t*_{2}, 24]), because the shortwave radiation is zero, the equation can be simplified as

The analytical solution of the land surface temperature at time *t* is

where *t*_{0} is the start time of the period.

In period 2 (*t* ∈ [*t*_{1}, *t*_{2}]), Eq. (7) is rewritten as

Multiplying *e*^{(B/A)t} on both sides, we have

We then integrate both sides at the interval [0, *t*]:

Because the shortwave solar radiation *Q*_{SW} is a piecewise function of time, the integration of Eq. (15) is divided into two parts: the first part is period 1 ([0, *t*_{1})) and the second part is from sunrise to time *t* ([*t*_{1}, *t*]); the current time *t* is the time between *t*_{1} and *t*_{2}. The integration of Eq. (15) is shown in appendix B. We derive the land surface temperature at time *t* (*t* ∈ [*t*_{1}, *t*_{2}]) as follows:

At time *t*_{1}, the results from Eqs. (12) and (16) are identical, which shows good agreement with the piecewise function. Last, the simple analytical solution of the daily temperature cycle is a combination of Eqs. (12) and (16).

The result is used as the thermal boundary condition for the csCFD model (Wang and Li 2016), which is modified on the basis of the traditional CFD model to simulate both the microscale environment and the city-scale environment. With csCFD, we can simulate the microscale environment and mesoscale environment simultaneously and obtain the city-scale environment (Wang and Li 2016; Hidalgo et al. 2008) with finer resolution. Because setting the sensible heat flux is consistent with the simple boundary conditions, the sensible heat flux derived from this temperature solution is implemented in the csCFD model (ANSYS 2011).

## 3. Data and results

### a. Evaluation against a benchmark experiment

The proposed sDCT boundary condition was first tested with benchmark field experimental data by Swaid and Hoffman (1989), who measured surface and air temperatures over a 9-month period. The meteorological conditions on 20 and 21 January 1987 and surface temperatures measured in bare soil were used for evaluation. The input data included solar radiation, air temperature, 10-m wind speed, and relative humidity (RH) provided by Swaid and Hoffman (1989, their Fig. 6). Shortwave radiation was calculated using the specific location and date at/on which the measurement took place; the measured radiation and the result from the analytical solution are shown in Fig. 2a. As explained by Swaid and Hoffman (1989), a shadow on the observation site after 1600 local solar time (LST) caused a significant reduction in solar radiation. Therefore, a compromise was needed because it is not possible to accurately predict both the time of sunset and sunrise with the proposed theoretical calculation. Since the reduction was just before sunset, we artificially set the shortwave radiation to zero from 1600 LST until sunset to conform to the actual conditions. This change will make sure that the calculated solar radiation approximates the actual radiation. Note that the calculated solar radiation was still somewhat overestimated after the peak hour, which leads to a slight overestimation of the calculated land temperature, as shown in Fig. 2b. This is because the ground surface to several centimeters below can be treated as a thermal mass in the calculation, which stores the energy that is overestimated during the calculation. At first, the overestimated energy was stored in the thermal mass and only slightly affected the air temperature. In the afternoon, the solar radiation decreased, but the energy in the thermal mass began to be released. This is also the reason that the highest air temperature often appears around 1400 LST rather than at noon. Because of the energy overestimation, the thermal mass releases more energy than in the observation, which leads to higher air temperatures in the calculation than in the observation.

Figure 2b shows the 24-h land temperature profile from the analytical solution and from the measurement on the selected day. The measured air temperature profile, which was used to compute the land surface temperature in the analytical solutions, is also shown for reference. The time at which the maximum temperature occurs was nearly the same in the calculation and in the observation. As mentioned previously, the calculated solar radiation in the simulation was higher than the observed radiation except at noon, which means that the land received more solar radiation in the simulation than in the observation. As a result, the time at which the surface temperature decreased in the afternoon differed slightly from that of the observation.

### b. Validating surface temperature from Hong Kong Observatory data

We also validated the analytical model by comparing the calculated temperature with the observed land surface temperature and other meteorological data obtained from the Hong Kong Observatory (HKO). All data with the exception of cloudiness were obtained from the King’s Park (KP) station, which is in an open area in a small park in the center of the city. The surface is soil. Data on cloudiness were obtained from the HKO station in Tsim Sha Tsui, because it is not recorded at the KP station. It is fortunate that the direct distance between the two stations is only ~1.2 km so that the cloudiness at HKO likely approximates that at the KP station, particularly on clear days. We picked a clear and sunny day for validation. The results from 8 March 2013 are shown here. The hourly average air temperature, RH, cloudiness, rainfall, solar radiation, and 10-m wind speed are used as the input data, and the hourly averaged grass surface temperature is used for validation. The albedo of the bare soil is used for the surface, and the emissivity is 0.9. The latitude of Hong Kong is 22.267°N.

Figure 3a shows the observed and calculated shortwave radiation. The measured shortwave radiation includes both diffuse and direct radiation. Figure 3b shows the observed hourly average RH and the wind speed, which is used to calculate the land surface temperature. The RH is low during the day and high at night, whereas the wind velocity is high during the day, particularly in the afternoon, and very low at night.

Figure 3c shows the calculated land temperature, the corresponding observed land temperature, and the observed air temperature profile at the KP station on 8 March 2013. The calculation starts from 7 March 2013 to avoid the influence of the initial conditions on the starting point. The calculation includes two days, and the result from the second day (8 March 2013) is used here for comparison. The calculated result agrees reasonably well with the observed grass temperature. The land temperature was nearly constant at night and increased rapidly after sunrise. The land temperature achieved its maximum at 1200 LST in both the calculated and observed results, when the shortwave solar radiation achieved its highest value. The surface temperature decreased during the afternoon as solar radiation decreased and less energy was transported to the surface. We also found that the agreement was much better at night than during the day. The calculation overpredicted the temperature, particularly at 0900 and 1600 LST. We believe that this overprediction is related to the calculation of the sensible heat flux. The method of calculating the sensible heat flux is simple; it is assumed to be proportional to the wind velocity.

### c. Different surface temperatures in urban and rural areas

Many factors contribute to the differences in surface temperature between urban and rural areas, including the thermal properties of ground materials, emissivity, albedo, wind velocity, and anthropogenic heat. Appendix C shows the contribution of most of these factors and the ability of the proposed method to predict it. By including anthropogenic heat, our model was able to predict the differences in the rural and urban surface heat fluxes. In the demonstration presented in the next section, the surface temperatures of the urban and rural areas were ultimately calculated with the average air temperature and wind velocity in the simulation. The adopted thermal properties and other parameters used for calculation are listed in Table 1. The anthropogenic heat in the rural area is zero, whereas the maximum anthropogenic heat was set to 100 W m^{−2} in the city. For the ideal city that we considered, we assumed that the anthropogenic heat varies with time in the urban area; this variability is the same as the diurnal variability of gross anthropogenic heat in spring in Tokyo, Japan, as shown in the study by Ichinose et al. (1999, their Fig. 6).

Figure 4a shows the air temperatures and calculated surface temperatures for the following case study. The hourly temperatures were calculated beginning at 0600 LST, before sunrise. It is assumed that the city is Hong Kong (22°15′0″N, 114°10′0″E) and the date is 8 March 2013. The initial air temperature is 24.85°C, and the initial surface temperature is 19.85°C (as in section 3b). On that day, the sun rose at 0639 LST and set at 1830 LST. The properties in Table 1 are used. The average wind velocity and air temperature at 10 m above the ground can also be obtained during the simulation. In our CFD simulation, the air temperature and wind velocity are not known at time step *n*. We used an explicit method in which the values at time step *n* − 1 are used to predict the surface temperature at *n*, which is one of the applications of the piecewise method. Although the rural air temperature is generally higher than that in urban areas during the day, when anthropogenic heat is included, the air temperature could be lower in rural areas than in the city (Yang and Li 2013). Such high temperatures as in Fig. 4a are also found in Higashiyama et al. (2016). Figure 4b shows the total heat flux released from the rural area and the urban surface, including the anthropogenic contribution. This heat flux is used as the boundary condition in the demonstration in the next section for calculation. As a result to anthropogenic heat, the heat flux released from the underlying surface of the city was always greater than the rural heat flux. The maximum heat flux in the city occurred at 1400 LST, as did the maximum surface temperature.

### d. Demonstration of its use in a city-scale CFD study

We consider a simple scenario in which the rural area is represented by a homogenous flat heating plate with different heat fluxes and there is no synoptic wind; the urban area is represented by a porous medium. Hang and Li (2010) and Hang et al. (2012) discussed the feasibility of treating the urban area as a porous medium for regular and irregular building clusters. Hu et al. (2012) used a porous medium and included heat in a simulation of the urban area. Wang et al. (2017) also discussed the possibility of using part of the city as a porous medium. It is assumed that the city is ideal and that its buildings are uniform cubes. The building density *λ*_{p} and frontal area density *λ*_{f} are both 0.25, which are characteristic of a city with low packing density, as defined by Grimmond and Oke (1999). It is further assumed that the city has an elevation of 50 m and that the city’s area is as shown in Fig. 5. A volume heat source was considered in the whole porous medium (Hu et al. 2012).

The flow is expected to typify urban heat island circulation (UHIC; Lu et al. 1997a,b), which typically forms during windless sunny days under two more preconditions: a stable PBL and a well-mixed convective PBL before the UHIC begins to form (Zeman and Lumley 1976). Therefore, we assume that the initial condition is stable and that the initial potential temperature lapse rate is 0.007 K m^{−1} before sunrise.

The computational domain is shown in Fig. 5 and is the same as in Wang and Li (2016). It is assumed that the city is infinitely wide in the spanwise direction and is homogenous and that the simulation domain can be simplified to two dimensions. The computational domain is 2950 m high and 110 km long. The city is located at the middle of the bottom of the domain with a length *D* of 10 km. The other bottom surface represents the rural area. The first grid in the *y* direction is 0.1 m, and the grid increases with a ratio of 1.05 until it reaches 50 m, which is the maximum grid spacing in the *y* direction. The *x* grid in the city area is 100 m, and the grid increases to its maximum of 500 m from the edge of the city to the rural area with a ratio of 1.05. The grid is the same as in Wang and Li (2016), which satisfies the requirements of the CFD guidelines for urban wind simulations (Franke 2006; Tominaga et al. 2008). The grid independence is shown in Wang and Li (2016) and is not represented here.

The boundary conditions are also shown in Fig. 5. The computational domain is symmetric with respect to the *y* axis or the line *x*/*D* = 0, and therefore to reduce the computation time only the right part of the domain is simulated. The hourly temperatures of the urban and rural surfaces are obtained with the new sDCT boundary conditions. The diurnal convective PBL develops after sunrise and dissipates slowly after sunset. The daytime heat flux in Fig. 4b is used for the diurnal convective PBL simulation shown in the next section. The fluxes are renewed every hour.

A third-order “quadratic upstream interpolation for convective kinematics” (QUICK) scheme is used for the convection terms, and second-order center differencing is used for the diffusion terms. The “pressure implicit with splitting of operator” (PISO) scheme is used for the pressure–velocity coupling method. The standard *k*–*ε* model (ANSYS 2011) is adapted for the turbulence model. The time step is 1 s; more details of the settings were described by Wang and Li (2016). The roughness height is not considered in the urban area because the city is represented as a porous medium, and the roughness of the rural area is calculated using the method described by Blocken et al. (2007). In the city area, the Darcy term, the Forchheimer term, and other related terms are used in the porous medium instead of the drag coefficient when considering the effects of the drag force.

### e. The development of urban heat island circulation

In the study of UHIC, the thermal environment and airflow before the final UHIC are also important. Figure 6 shows the bihourly development of the UHIC with both vector maps and contour maps. Figure 7 shows the hourly averaged horizontal velocity profile along the city center (*x*/*D* = 0), the city edge (*x*/*D* = 0.5), and rural areas (*x*/*D* = 1), respectively. When these two figures are examined at the same time, the evolution of the UHIC becomes clear.

In the beginning, the circulation is small; it appears only near the city edge and is nearly zero in the middle. The horizontal temperature gradient is then replaced by the advection of colder air due to the development of winds toward the city center, which in turn creates perturbations in the temperature field. The velocity increases from the city center to the city edge. A pressure gradient is finally created and enlarges the circulation cell. As a result, the horizontal velocity is nearly zero in the city center and the vertical velocity is also nearly zero at the city edge during the first few hours, which drives the wind that develops toward the city center (Figs. 7a,b). As the circulation cell gets larger, the center of the cell moves toward the *y* axis and leads to some downward vertical velocity in front of the cell before the cell center reaches the *y* axis, which explains why the vertical velocity is negative at 1000 LST (Fig. 7b). For the same reason, the horizontal velocity at the city edge is greater than that in the rural areas at 1000 LST (Fig. 7a). Meanwhile, the air rises above the city because of the urban geometry.

After 1000 LST, the air converges in the lower part of the urban boundary layer near the city center; it then rises above the city center and diverges at the upper part of the PBL, thus forming a complete air circulation driven by the effective sensible heat flux. This circulation cell is the UHIC. This result agrees with the development of the circulation given by Delage and Taylor (1970).

From 1200 to 1800 LST, the effective heat flux increases as a result of the greater thermal storage capacity and anthropogenic heat of the city. The vertical and horizontal velocities also increase, which demonstrates the increase in the strength and size of the UHIC (Figs. 6c–e). The difference in the horizontal velocities between the city edge and the rural area is small in the afternoon (Fig. 7a), after the UHIC has formed. The maximum height of the vertical velocity also increases during the afternoon.

Figure 7c shows the hourly average potential temperature profile along the city center (*x*/*D* = 0), the city’s edge (*x*/*D* = 0.5), and the rural area (*x*/*D* = 1) during the day. The thermally unstable surface layer and the overshoot area recover over time. (The same phenomenon can be found in Fig. 6.) The UHIC continues to rise to the base of the inversion layer, leaving a convective layer below. The convective PBL develops more quickly in the city than in the rural area, and the overshoot height is higher in the city than in the rural area. The height decreases from the city center to the rural area. In the convective PBL, the potential temperature in the city is always higher than that in the rural area.

The height of the mixing layer (the mixing height) is one of the key parameters for the scaling and convective PBL, but the definition and calculation of the mixing height are not uniform in all studies (Stull 1988; Seibert et al. 2000; Garratt 1994). Seibert et al. (2000) noted that the definition of the mixing height differs for the convective PBL and for a neutral or stable boundary layer. In the convective urban boundary layer, “the height where the heat flux gradient reverses its sign” is the value closest to the thermodynamic definition of the convective urban boundary layer mixing height (Seibert et al. 2000). This method is widely used for measurement at high temporal frequencies when the turbulence characteristics are available. For situations in which the turbulent fluxes are unavailable because of limitations in the sampling frequency, Heffter (1980) provided a method to estimate the mixing height from profiles of potential temperature. This method has been widely used as the mixing height and is easily obtained, and it is used here. In this study, the base of the inversion layer must fulfill two requirements: 1) the potential temperature lapse rate must be no less than 0.005 K m^{−1} and 2) the potential temperature at the top of the inversion layer should be at least 2 K higher than that at the base of the inversion layer. The mixing height is the height above the base of the inversion layer with a potential temperature that is 2 K higher than that at the base of the inversion layer.

The predicted mixing heights of the city center, city edge, and rural area at each hour are determined from the hourly average potential temperature profiles. The mixing height at the city center is higher than that at the city edge and in the rural area during the day because more heat is released from the city than from rural areas. The highest mixing height is observed in the late afternoon (1800 LST). When the porous city is not considered, the highest mixing height is observed at 1600 LST. These findings agree well with measured values (Dupont et al. 1999).

Understanding the evolution of the UHIC is very important in the urban environment, especially when a serious pollution episode occurs, such as the haze that is frequently seen in some large cities across China. Such serious haze episodes often occur when the background wind is calm, that is, when the UHIC is present. Studies have mainly focused on situations in which the UHIC has formed, referred to as a quasi-steady situation, but not on the UHIC’s evolution. A better boundary condition that includes a daily cycle could improve our understanding of urban airflows in these situations and help us to identify ways to eliminate them. The simulation also shows the development of the thermal and wind environments, not just in the urban canopy, but also in the PBL. Meanwhile, because the simulation domain is large and a fine grid resolution is required, the use of a simple surface boundary condition for city-scale environmental analysis is more effective than the use of a coupled model.

## 4. Discussion

There are several methods for estimation of the surface temperature. Some focus on the entire energy budget of the whole city or a larger area (e.g., Masson 2000), and some focus on the microscale temporal distribution of the surface temperature in cities or building clusters (Krayenhoff and Voogt 2007; Yang and Li 2015). They all must numerically solve the partial differential equations of the energy balance with iteration. We attempted to solve the equations analytically to save computation time during the simulation, especially for the fine grid, when the time step is small. One advantage of our proposed method is its efficiency as a fast method with a semianalytical approach. In theory, if we assume that the conductive heat transfer through all surfaces in an urban area is one-dimensional, we can also integrate it with a 3D radiation model and obtain the temperature at each surface point in the city. However, storage difficulties still prevent the existing 3D radiation method from simultaneous simulation of thousands of buildings in a city.

In theory, the ground temperature could also be predicted by simplifying the urban canopy models with some modifications. For example, Kanda et al. (2005) and Kawai et al. (2009) could predict surface temperatures at each face of the urban canopy (i.e., roof, floor, and four vertical walls) well for the city with iteration. The method proposed here could also predict the surface temperatures of all surfaces of the building explicitly, and the process is much simpler: by changing the property of the underlying surface, the convective heat transfer coefficient, the albedo, and other parameters. When we need to consider a real city—for example, a large amount of buildings or the heterogeneity of the city (at present we assumed that the buildings are a regular array as in the urban canopy model), this method could calculate the surface temperature much faster than can the urban canopy models and other models.

Williamson (2004) also derived an analytical solution, but this derivation linearized the longwave radiation to a first-order formula, *B*(*T*_{s} − *T*_{a}), in which the gradient of the linearized function *B* is a constant. However, *B* is actually related to the air temperature, which is always changing. Meanwhile, the sky temperature instead of the air temperature should be considered when calculating the longwave radiation. Silva et al. (2009) indicated that the sky temperature could be calculated from the dry air temperature and the RH. The formula for longwave radiation is improved here by including the sky temperature and improving the gradient term. After this improvement, the calculated surface temperature is improved, especially at night.

The method is currently only validated for sunny days, and its feasibility for rainy days requires further validation. Meanwhile, the heat transfer coefficient of the sensible heat flux is now simply approximated. We are aware that there are many choices of the heat transfer coefficient; at some point, a sensitivity study should be performed to choose the best heat transfer coefficient. Another issue at present is that the method can consider only the ground surface and roofs of buildings. For the sidewalls of buildings, the radiation terms require modification to include the effects of shadows and reflections that result from the different height–width ratios.

## 5. Conclusions

We have proposed an sDCT boundary condition that provides a 24-h temperature change to consider the main factors of the energy budget balance for surfaces, including solar radiation, longwave radiation, sensible heat flux, soil heat flux, and anthropogenic heat. The advantages of this boundary condition include the use of temperature as a simplified analytical solution for the energy budget and the fact that computational time is reduced because there is no need to couple an energy balance model with the traditional CFD model.

We evaluated the performance of the sDCT boundary condition with a benchmark experiment from the literature and with data from the Hong Kong Observatory. The simulated results show good agreement with the observed surface temperatures. The contributions of the thermal properties of materials, albedo, the sky-view factor, and anthropogenic heat are also discussed. Our model closely predicts that large thermal storage delays the occurrence of the maximum surface temperature. The new sDCT is evaluated for calculating the daily circulation that is induced by the urban heat island by using csCFD when the synoptic wind is zero and the city is 10 km wide.

## Acknowledgments

This work was supported financially by a Research Grants Council Collaborative Research Fund project grant (HKU9/CRF/12G) from the Hong Kong SAR government. We thank Mr. Lei Hao for his help in obtaining the analytical solution of the equations. This work was also supported by the Natural Science Foundation of Guangdong Province with Grant 2016A050503035. The Hong Kong weather data were obtained from the Hong Kong Observatory. All other data and the computer code used in this paper can be obtained from Xiaoxue Wang.

### APPENDIX A

#### Description of Each Term in the Surface Energy Balance Equation

The *net radiation Q*_{NR} is the sum of all shortwave and longwave radiation, which includes the incoming and outgoing shortwave radiation and incoming and outgoing longwave radiation:

*Shortwave radiation* has been discussed in section 2.

*Longwave radiation* includes that emitted from the surface and that emitted from the sky . Both are dependent on the emissivity:

where *T*_{sky} (K) is the sky temperature, *T*_{l} (K) is the land temperature, *T*_{a} (K) is the air temperature, *F*_{sky} is the sky-view factor (SVF) of the area of interest, *E* is the emissivity, and *σ* is the Stefan–Boltzmann constant (5.67 × 10^{−8} W m^{−2} K^{−4}).

The longwave radiation can thus be calculated as

The sky temperature is related to the RH and the dewpoint. Silva et al. (2009) showed a method to calculate the sky temperature from the dry-bulb air temperature *T*_{a} and the RH:

the dewpoint is calculated from

where

*a* = 17.27, and *b* = 237.7. The temperature is always expressed in kelvins.

To obtain an analytical solution, Eq. (A4) is linearized by making the appropriate approximation. First, we assume that *E*_{l} ≈ *E*_{a} so that Eq. (A4) is rewritten as

where

is the gradient of the linearized function and *T*_{sky} is calculated by *T*_{a}(*e*_{at})^{0.25}.

The gradient of the linearized function can be rewritten as

Note that *h*_{a} is not sensitive to the air temperature. When the air temperature is between −35° and 65°C, and the difference between the air temperature and the land surface temperature is less than 50°C, the deviation of *h*_{a} is less than 0.25, as shown in Fig. A1. The gradient of the linearized function and the longwave radiation can thus be simplified as

where *e*_{at} = 0.004*T*_{dew} + 0.8.

The sensible heat flux is given by

where *U*_{10} is the wind velocity at a reference height (in this case, 10 m above the ground), *ρ*_{a} = 1.225 kg m^{−3} is the air density, *C*_{p} is the specific heat capacity of air, and *C*_{H} is the transfer coefficient for heat given by

with *κ* being the von Kármán constant, *z*_{r} being the reference height at which *U* is measured, and *z*_{0} being the roughness length (here 0.1 m).

The heat transfer coefficient can also be calculated from the experimental equations. One of the experimental equations that is based on the “surface thermal time constant” model (Swaid and Hoffman 1989) is

The form of the heat transfer coefficient is not unique. Many other kinds of equations can be used to calculate the sensible heat flux. We have chosen this method because it is easier to implement, particularly for analytical solutions. It will be improved by more accurate calculations in the future.

*Soil heat flux* depends on changes in the soil’s temperature and heat capacity:

where *C*_{l} represents the volumetric soil heat capacity and ∂*T*/∂*t* is the rate of soil temperature change. The volumetric soil heat capacity (J m^{−3} K^{−1}) can be expressed as

in which *ρ*_{l} is the sandy soil density (1.6 × 10^{3} kg m^{−3}) and *c*_{l} is the specific heat of the soil (0.35 × 10^{3} J kg^{−1}K^{−1}).

By considering a thin layer at the surface with a depth scale *h*_{l}, we have

The choice of the depth scale must be made very carefully. Stull (1988) proposed a method to determine *h*_{l}:

where *P* = 1 day (in seconds), *ν*_{g} = *k*_{g}/*C*_{l} is the soil thermal diffusivity, and *k*_{g} is the molecular conductivity of the soil. The value of thermal diffusivity is shown in appendix C-8 in Stull (1988). Here, the soil thermal diffusivity is 0.38 × 10^{−6} m^{2} s^{−1}.

When the urban canopy layer is not considered, the total anthropogenic heat *Q*_{AN} may be considered in the surface heat balance. The literature contains considerable data on anthropogenic heat (e.g., Ichinose et al. 1999).

### APPENDIX B

#### The Integration Process of Eq. (14)

When *t* ∈ [*t*_{1}, *t*_{2}], Eq. (14) can be written as

Before sunrise, *Q*_{s} is zero. The surface temperature at sunrise *t*_{1} can be obtained by integrating Eq. (14) from 0 to the time of sunrise. At any time between 0 and *t*_{1}, the cosine function of time is the same as the cosine function of *t*_{1} ( = constant). Using time *t*_{1} to represent the cosine function in Eq. (14) and integrating the left-hand side of Eq. (B1), we have

Integration of the last term in Eq. (B3) yields

The equation can be rewritten as

### APPENDIX C

#### Factors that Contribute to the Surface Temperature Difference

To determine the contribution of these factors, the same air temperature and RH observed on 8 March 2013 at the KP station were used as the input, and the value of each factor was changed separately for our study. The wind velocity was also assumed to be constant so that it would not influence the temperatures of different surfaces or areas. The land temperatures from 0000 LST on 7 March 2013 to 0000 LST on 9 March 2013 were calculated. To avoid the influence of the initial conditions, only the results from the second day are shown in the following results.

First, urban and rural areas have different underlying surface materials. Rural areas are often covered in grass or bare soil, whereas cities are usually covered with concrete. The different materials have their own thermal properties, including density, specific heat, and thermal diffusivity, and ultimately cause different surface temperatures. Table C1 shows the thermal physical properties of four common materials in typical urban and rural areas.

Figure C1a shows the 24-h temperature of surfaces made of four different materials. The SVF for each case is assumed to be 1.0, and city structures are not considered. The temperature difference among the four materials was very small at night but increased during the day. Because of the greater thermal storage of concrete, the maximum temperature of the concrete surfaces occurs later than that of other surfaces, and the temperature is higher than those of the other surfaces at night. The differences in the thermal properties between sandy soil and clay were small. The differences in temperature were a result of the combination of thermal storage and emissivity. Among the soil surfaces, the temperature variation of the peat soil was the greatest, whereas that of the sandy soil was the least.

Second, the different arrangements of a city’s buildings lead to differences in albedo and SVF that influence the solar radiation at the surface (Oke 1978; Yang and Li 2013). Yang and Li (2013) showed that the arrangement of buildings could influence both the albedo and the SVF of the road, and calculated those values under different typical plan area indices, as shown in Table C2. To evaluate the contribution of the albedo and SVF, we calculated the surface temperature for the day used above, 8 March 2013, with a concrete surface material and a constant wind velocity, resulting in only albedo and SVF changes.

Figure C1b shows the daily surface temperature variations when the albedo and SVF in a city change as a result of different plan area indices. When the plan area index is 0, the city has no buildings, and therefore the SVF is 1.0 and the initial albedo is 0.41, as given by Yang and Li (2015). When the arrangement of buildings is denser, the albedo and the SVF both decrease, which means that less radiation heat is released and the surface temperature increases during the day. The warming effect achieves its peak when the plan area index is around 0.65. From this point, when the buildings are closer, the roof behaves like a surface and reflects much of the radiation, and the albedo begins to increase. Consequently, although the SVF decreases, the maximum surface temperature also decreases.

## REFERENCES

*Fluent 14.0 User’s Guide*. ANSYS, Inc., 2498 pp.

*Air Pollution Meteorology and Dispersion*. Oxford University Press, 310 pp.

*Fourth Int. Symp. on Computational Wind Engineering*, Yokohama, Japan, Int. Association for Wind Engineering, 529–532, http://www.iawe.org/Proceedings/CWE2006/TA4-01.pdf.

*The Atmospheric Boundary Layer*. Cambridge University Press, 316 pp.

*Boundary Layer Climates*. Methuen and Co., 435 pp.

*An Introduction to Boundary Layer Meteorology*. Kluwer Academic, 670 pp.

*Seventh WRF Users’ Workshop*, Boulder, CO, UCAR, 5.6, http://www2.mmm.ucar.edu/wrf/users/workshops/WS2006/abstracts/Session05/5_6_Tewari.pdf.

## Footnotes

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).