The Microscale Urban Surface Energy (MUSE) Model for Real Urban Application

: Urban atmospheric environmental issues are commonly associated with the physical processes of urban surfaces. Much progress has been made on the building-resolving microscale atmospheric models, but a realistic representation of the physical processes of urban surfaces on those models is still lacking. This study presents a new microscale urban surface energy (MUSE) model for real urban meteorological and environmental applications that is capable of representing the urban radiative, convective, and conductive energy transfer processes along with their interactions, and that is directly compatible with the Cartesian grid microscale atmospheric models. The physical processes of shadow casting and radiative transfers were validated on an analytical accuracy level. The full capability of the model in simulating the three-dimensional surface heterogeneities in a real urban environment was tested for a hot summer day in August 2016 using the ﬁeld measurements obtained from the Kongju National University campus, South Korea. The validation against the measurements showed that the model is capable of predicting surface temperatures and energy balance ﬂuxes in a patch scale at the heterogeneous urban surfaces by virtue of the interactive representation of the urban physical processes. The excellent performance and ﬂexible grid design emphasize the potential capabilities of the MUSE model for use in urban meteorological and environmental applications through the building-resolving microscale atmospheric models, such as computational ﬂuid dynamics (CFD) and large-eddy simulations (LES).


Introduction
Microscale atmospheric models, such as computational fluid dynamics (CFD) and large-eddy simulations (LES), have been widely used to simulate urban flow, turbulence, and pollutant dispersion in urban environments (e.g., [1][2][3][4][5][6]). Many previous studies focused on understanding urban flow behaviors and associated dispersion characteristics using extremely simplified urban geometries. For example, Kim and Baik [7] investigated the flow and dispersion characteristics in a single canyon by changing the canyon aspect ratio (that is, the ratio of building height to canyon width). Lin et al. [8] and Ramponi et al. [9] examined the influence of street width on local ventilation, and Hang et al. [10] demonstrated a significant impact of building height change on ground-level pollutant concentrations using simple cubic array type geometries.
Recently, the models' applications have rapidly extended to answer urban atmospheric environmental issues of real urban environments such as ventilation for wind and thermal comfort (e.g., [11][12][13][14][15]). The change is attributed to urban geographical information data and computing resources [16]. The geographic information system (GIS) data that characterize the morphological features of urban surfaces in many cities for buildings, streets, trees, and others have been readily available (e.g., [17]). Moreover, the development of high-performance computers has made CFD and processes and comparisons against the field measurements. A summary and conclusions follow in Section 5.

Grid Representation of Urban Surfaces
Real urban surfaces consist of various artificial and natural obstacles such as buildings, roads, trees, and grasses, and their different combinations create unique urban geometries. The realistic representation of the urban surfaces is a prerequisite in parameterizing urban physical processes for microscale flow and dispersion simulations. The MUSE model represents the urban surfaces with three-dimensional structures of grid cells on a Cartesian grid system (Figure 1a). Each grid cell has six "patches" of two horizontal (top and bottom) surfaces and four vertical (east, west, south, and north) wall surfaces. Among them, the patches facing the atmosphere are defined as active patches, on which surface energy exchange processes are calculated between the patch and the corresponding atmospheric model grid. and are the air temperature and wind speed at the adjacent atmospheric grid cells, respectively.
The MUSE model is designed to the three-dimensional grid structure that is compatible directly with general microscale atmospheric models (Figure 1b). Thus, the interaction between microscale atmospheric models and the MUSE model can be calculated efficiently with both "online" or "offline" modes. Table 1 lists the meteorological forcing variables and the predicted variables of the MUSE model. The online computation means that the time integration of the MUSE model and any microscale atmospheric model can be undertaken in a simultaneous step with each other. In other words, the MUSE model predicts the surface energy balance fluxes and surface temperatures at each patch given meteorological fields by the microscale atmospheric model, and the microscale atmospheric model predicts three-dimensional meteorological fields using the surface energy balance fluxes or surface temperatures as a bottom boundary condition. In offline computation mode, the predictions of the MUSE model and the microscale atmospheric model are conducted sequentially no matter which model runs first. When meteorological forcing measurements are available, it is recommended to use the offline mode computation in reproducing the thermal heterogeneities of the entire urban area. Generally, current microscale atmospheric models cannot provide radiative and thermal forcing of the MUSE model as accurate as measurements due to lack of corresponding physical representation. When it comes to dynamic forcing (that is, wind fields), measurements have limitations in representing spatial heterogeneities of complicated urban surfaces due to low The MUSE model is designed to the three-dimensional grid structure that is compatible directly with general microscale atmospheric models ( Figure 1b). Thus, the interaction between microscale atmospheric models and the MUSE model can be calculated efficiently with both "online" or "offline" modes. Table 1 lists the meteorological forcing variables and the predicted variables of the MUSE model. The online computation means that the time integration of the MUSE model and any microscale atmospheric model can be undertaken in a simultaneous step with each other. In other words, the MUSE model predicts the surface energy balance fluxes and surface temperatures at each patch given meteorological fields by the microscale atmospheric model, and the microscale atmospheric model predicts three-dimensional meteorological fields using the surface energy balance fluxes or surface temperatures as a bottom boundary condition. In offline computation mode, the predictions of the MUSE model and the microscale atmospheric model are conducted sequentially no matter which model runs first. When meteorological forcing measurements are available, it is recommended to use the offline mode computation in reproducing the thermal heterogeneities of the entire urban area. Generally, current microscale atmospheric models cannot provide radiative and thermal forcing of the MUSE model as accurate as measurements due to lack of corresponding physical representation. When it comes to dynamic forcing (that is, wind fields), measurements have limitations in representing spatial heterogeneities of complicated urban surfaces due to low instrumentation Atmosphere 2020, 11, 1347 4 of 26 density. The neutral simulation of microscale atmospheric models might be a better choice because of the models' capability in resolving surface drags of buildings. Overall, the colocated grid design makes the MUSE model flexible to meteorological forcing in any combination with microscale atmospheric models and measurements. Various physical processes associated with the complicated urban surfaces interact, resulting in the formation of unique urban microclimates (e.g., [48][49][50][51][52][53][54][55][56]). The MUSE model includes the urban physical processes that are shortwave and longwave radiative transfers within the canyon, turbulent sensible heat exchanges at each patch surface, and conductive heat transfer within multiple interior subsurface layers over an urban area ( Figure 2). It also includes a shadow submodel and the formulations of three-dimensional patch-patch view factor estimation for the radiative transfer parameterization. All of those physical processes interact simultaneously among them during the model integration, leading to realistic reproduction of the urban surface heterogeneities in surface energy balance fluxes and surface temperature.
The urban surfaces of buildings and roads are made with diverse construction materials such as concrete, brick, and asphalt. Thus, the urban surface simulation needs to characterize the urban surfaces realistically in association with urban physical processes. The MUSE model characterizes urban surfaces using the physical property parameters assigned at each patch of the simulation domain. Table 2 lists the physical property parameters of the MUSE model used to represent the physical properties of urban surfaces. The surface roughness lengths for momentum and heat are the aerodynamic property parameters associated with turbulent momentum and heat exchanges, respectively. The surface albedo and emissivity are the radiative property parameters of urban surfaces, and the thermal conductivity and volumetric heat capacity are the thermal property parameters of the interior subsurfaces. Each patch has all the physical property parameters assigned according to urban surface materials.
Atmosphere 2020, 11, 1347 5 of 26 Atmosphere 2020, 11, x FOR PEER REVIEW 5 of 26 Figure 2. Configuration of the subsurface layers at each patch and associated physical processes in the MUSE model. , denotes the subsurface temperature at the layer of the patch , and , is the conductive heat flux at the vertical level of the patch . ↓ and ↓ denote the downward direct and diffuse shortwave radiation, respectively, and ↓ is the downward longwave radiation. and are the air temperature and wind speed at the adjacent atmospheric grid cells, respectively.
The urban surfaces of buildings and roads are made with diverse construction materials such as concrete, brick, and asphalt. Thus, the urban surface simulation needs to characterize the urban surfaces realistically in association with urban physical processes. The MUSE model characterizes urban surfaces using the physical property parameters assigned at each patch of the simulation domain. Table 2 lists the physical property parameters of the MUSE model used to represent the physical properties of urban surfaces. The surface roughness lengths for momentum and heat are the aerodynamic property parameters associated with turbulent momentum and heat exchanges, respectively. The surface albedo and emissivity are the radiative property parameters of urban surfaces, and the thermal conductivity and volumetric heat capacity are the thermal property parameters of the interior subsurfaces. Each patch has all the physical property parameters assigned according to urban surface materials.  Configuration of the subsurface layers at each patch and associated physical processes in the MUSE model. T k,I denotes the subsurface temperature at the layer k of the patch I, and F k,I is the conductive heat flux at the vertical level k of the patch I. S dir↓ and S di f ↓ denote the downward direct and diffuse shortwave radiation, respectively, and L di f ↓ is the downward longwave radiation. T a and U a are the air temperature and wind speed at the adjacent atmospheric grid cells, respectively.

Shadow Model
A shadow model is developed as an ingredient of the MUSE model based on a geometric ray casting approach. In real urban environments, the shadow geometry of each building depends on the relative position of the surrounding buildings and the sun. The shadow model at first step computes the shadow lengths at each rooftop patch of buildings cast on horizontal and vertical surfaces, and then it determines the sunlit/shaded condition at every patch ( Figure 3). The shadow length by the patch I of a building rooftop cast on a horizontal surface (l hor,I ) can be analytically determined by Atmosphere 2020, 11, 1347 where h I is the height of the patch I of a building rooftop and θ z is the solar zenith angle. The vertical shadow length at a certain point is computed along the horizontal shadow path (l vert,I ) by l vert,I = l hor,I − l P,I tan where l P,I (= ∆x 2 + ∆y 2 ) is the horizontal distance between the point and the patch I. Here, the point location is defined by a relative distance separated from the patch I with ∆x and ∆y. Therefore, the shaded patches by the patch I of a building are determined as the patches that exist on the horizontal and vertical shadow paths. The MUSE model inspects every patch on the horizontal shadow path moving forward from the location of patch I with a small distance step of (dx, dy) = (L cos θ a , L sin θ a ). Here, θ a is the solar azimuth angle and L is the step size of length.
A shadow model is developed as an ingredient of the MUSE model based on a geometric ray casting approach. In real urban environments, the shadow geometry of each building depends on the relative position of the surrounding buildings and the sun. The shadow model at first step computes the shadow lengths at each rooftop patch of buildings cast on horizontal and vertical surfaces, and then it determines the sunlit/shaded condition at every patch ( Figure 3). The shadow length by the patch of a building rooftop cast on a horizontal surface ( , ) can be analytically determined by where ℎ is the height of the patch of a building rooftop and is the solar zenith angle. The vertical shadow length at a certain point is computed along the horizontal shadow path ( , ) by where , (= ∆ + ∆ ) is the horizontal distance between the point and the patch . Here, the point location is defined by a relative distance separated from the patch with ∆ and ∆ . Therefore, the shaded patches by the patch of a building are determined as the patches that exist on the horizontal and vertical shadow paths. The MUSE model inspects every patch on the horizontal shadow path moving forward from the location of patch with a small distance step of ( , )=( cos , sin ). Here, is the solar azimuth angle and is the step size of length.
where is the latitude of the model domain, is the solar declination angle, and is the solar hour angle. The value of is set to increase clockwise from the due north. The solar zenith angle (θ z ) and solar azimuth angle (θ a ) at the center of model domain are calculated by cos θ z = sin φ sin δ + cos φ cos δ cos ω t ,

Patch View Factors
where φ is the latitude of the model domain, δ is the solar declination angle, and ω t is the solar hour angle. The value of θ a is set to increase clockwise from the due north.

Patch View Factors
The three-dimensional view factors are necessary for calculating radiative energy exchange processes within urban canyons. Each patch sees other patches and/or the sky. When considered for two patches I and J, the patch-patch view factor (ψ I→J ) denotes the fraction of the radiation reached at the patch I to that emitted from the patch J [57]. Similarly, the patch-sky view factor (ψ I→sky ) at the patch I indicates the fraction of radiation reached from the sky. Generally, analytical methods can be used for calculating ψ I→J in simple urban geometries (e.g., [58][59][60]). However, they do not apply to complicated morphologies in real urban environments. Some previous studies have used ray tracing algorithms (e.g., [40,42,43]). When a large number of photons tracks in the algorithm, it can calculate the view factors with high computation accuracy [61,62]. However, it inherently requires tremendous computing time to compute the view factors over hundreds of buildings for real urban simulations. Lee et al. [63] suggested a new numerical computation method for patch-patch view factor calculation to solve the limitation. The method was derived from the analytical view factor formulation, which can be applied to arbitrarily facing patches and is not limited to any specific patch orientation.
In this study, the analytically-based numerical method by Lee et al. [63] is used for perpendicular and parallel patches to calculate the patch-patch view factors and the patch-sky view factors ( Figure 4). When the two patches I and J are perpendicular to each other with a distance of R, the patch-patch view factor (ψ I→J ) can be computed as follows: where ψ E⊥W I→J and ψ S⊥N I→J are the patch-patch view factors from the horizontal patch I to the vertical patch J on east/west and north/south walls, respectively. ∆X = x I − x J , ∆Y = y I − y J , and ∆Z = z I − z J are the directional distance between the two patches (R = . ∆A J is the area of patch J (∆A J = ∆y J ∆z J for the east or west wall, ∆A J = ∆x J ∆z J for the north and south wall, ∆A J = ∆x J ∆y J for horizontal surface). When the two patches I and J are parallel to each other with a distance of R, the patch-patch view factor (ψ I→J ) can be computed as follows: where ψ E//W I→J and ψ S//N I→J are the patch-patch view factors from the vertical patch I to the vertical patch J on east/west and north/south walls, respectively. From simple algebra relation of view factor, the patch-sky view factor (ψ I→sky ) at a horizontal or vertical patch can be calculated as follows: where N denotes the number of patches that are facing the patch I. Note that N is determined for vertical patches by the patches that are above a horizontal plane from the patch I.
the patch-sky view factor ( → ) at a horizontal or vertical patch can be calculated as follows: where N denotes the number of patches that are facing the patch . Note that N is determined for vertical patches by the patches that are above a horizontal plane from the patch . wall-wall patches (modified from Lee et al. [63]). ∆A I and ∆A J are the areas of the patches I and J, respectively. ∆X, ∆Y, and ∆Z are the distances between the two patches in each directions, and R is the distance between the two patches. n I and n J are the unit vectors perpendicular to the patches I and J, respectively.

Shortwave Radiation
The shortwave radiative fluxes reaching every patch are composed of the direct and diffuse shortwave radiation. The direct shortwave radiative flux at the patch I (S dir I ) is calculated by where S dir↓ is the direct shortwave radiation at the horizontal surface from the sky, and c I is the solar incidence angle on the patch I, and f sunlit is the surface sunlit fraction. f sunlit is determined by the shadow model as a value of 1 (sunlit) or 0 (shaded) at each patch. The solar incidence angle on the patch I (c I ) is calculated as where β I and γ I are the slope angle and azimuth angle of the patch I, respectively. The value of γ is set to increase clockwise from the due north as θ a . This equation can apply for arbitrarily tilted surfaces. In the MUSE model, β has 0 for horizontal (roofs and roads) patches or π 2 for vertical (walls) patches. Thus, the resultant S dir I is identical to the downward shortwave radiation from the sky (S dir↓ ) for horizontal sunlit patches and that is S dir↓ cos c I cos θ z for vertical sunlit patches. It is assigned to 0 W m −2 for the other shaded patches. On the other hand, the diffuse shortwave radiative flux at the patch I under the assumption that the downward shortwave radiation from the sky (S di f ↓ ) is isotropic. The reflected shortwave radiative flux reaching the patch I from the other patch J can be expressed as where S R I,m+1 is the reflected shortwave radiative flux by the patch I reaching from the other patches J at (m + 1)th reflection, α I is the surface albedo of the patch I, N is the number of the patches J interacting with the patch I, and S R J,m is the shortwave radiative flux reflected from the patch J after mth reflection. Here, the Lambertian reflection is assumed at every surface. The initial reflected shortwave radiative where M is the number of reflection considered for radiative energy trapping within urban canyons, and three times reflection applies to the computation of shortwave radiation in the MUSE model. Meanwhile, the reflection computation is a time-and memory-consuming process. The MUSE model applies the compressed row storage scheme suggested by Yang and Li [42]. The number of patches interacting with the patch I limits with a small fraction of the entire patches, meaning that most of ψ I→J is zero. The compressed row storage scheme takes only nonzero values into account, which results in a significant reduction of computational cost for calculating multiple reflection. Because the scheme is simply to rearrange the nonzero view factor patches for computational efficiency, it causes no additional computational error. The MUSE model assumed a minimum threshold value of 1.2 × 10 −20 to further enhance computational efficiency, which leads to a negligible difference in surface temperature.

Longwave Radiation
Unlike shortwave radiation, the MUSE model considers the longwave radiation emitted from the whole patches as well as the downward longwave radiation from the sky. The diffuse longwave radiative flux reaching the patch where L di f ↓ is the downward longwave radiative flux from the sky. The longwave radiation emitted from the patch I (L I ) is calculated following the Stefan-Boltzmann equation as follows: where ε I and T I are the surface emissivity and surface temperature of the patch I, respectively, and σ is the Stefan-Boltzmann constant (=5.67 × 10 −8 W m −2 K −4 ). The computation of the reflection process of longwave radiation is similar in shortwave radiation except that the longwave radiation emitted at each patch is considered along with the diffuse longwave radiation from the sky. The reflected longwave radiative flux reaching the patch I from the other patch J can be expressed as where L R I,m+1 is the reflected longwave radiative flux from the patch I reaching from the other patches J at reflection (m + 1)th reflection, ε I is the surface emissivity of the patch I, N is the number of patches J, and L R J,m is the longwave radiative flux reflected from the patch J after (m)th reflection. The initial reflected longwave radiative flux at each patch can be calculated by L R I,1 = (1 − ε I ) L di f I + L I . The resultant net longwave radiative flux absorbed at the patch I (L * I ) can be expressed as where M is the number of reflection considered for radiative energy trapping within urban canyons, and one reflection applies to the computation of longwave radiation in the MUSE model.

Horizontal Roof and Road Patches
The MUSE model applies the Monin-Obukhov similarity relation for calculating turbulent sensible heat fluxes at the horizontal roof and road patches. Even though it is difficult to justify the applicability of Monin-Obukhov similarity relation on a theoretical basis, many mesoscale and microscale modeling studies apply the formulation with a practical purpose to determine physical fluxes at urban surfaces (e.g., [64,65]). The turbulent sensible heat flux at the horizontal patch I (H I ) is given following Lee and Park [64] by where ρ is the air density (kg m −3 ), c p is the specific heat capacity at constant pressure (J kg −1 K −1 ), u * is the friction velocity (m s −1 ), and T * is the temperature scale (K). κ is the von Kármán constant, z 0,I and z 0T,I are the surface roughness lengths for momentum and heat, respectively, Ri B is the bulk Richardson number (= gz∆T θU a 2 ), and F h is the atmospheric stability correction function. T I is the surface temperature of the patch I, and T a and U a are the temperature and wind speed at the adjacent atmospheric grid of the patch I, respectively.

Vertical Wall Patches
The empirical formulation based on Rowley et al. [66] has been widely applied for calculating the convective heat flux at vertical surfaces (e.g., [40,67]) because the use of the Monin-Obukhov similarity relation for vertical surface fluxes calculation is arguable. In the MUSE model, the turbulent sensible heat flux at the vertical patch I (H I ) is given by where T I is the surface temperature of the vertical patch I, and T a is the air temperature at the adjacent atmospheric grid cell of the patch I. The convective heat transfer coefficient (h c ) is given by where U a = √ U 2 + V 2 + W 2 is the wind speed at the adjacent atmospheric grid cell of the patch I. The empirical formulation was obtained from in situ measurements [66].

Subsurface Heat Conduction
Every active patch consists of multiple subsurface layers to represent conductive heat transfers in buildings and grounds. This process is one of the essential physical processes in the formation of urban heat islands [68]. The MUSE model assumes that the conductive heat process occurs in a one-dimensional column ( Figure 2). From the conservation of energy, the media temperature at each subsurface layer of the patch I can be calculated with an explicit finite differencing by where T n k,I is the media temperature of the subsurface layer k of the patch I at nth time step, C k,I is the volumetric heat capacity of the subsurface material, ∆z k is the depth of each sublayer, and ∆t is the integration time interval. F n k+1,I is the conductive heat flux at the interfacing level between the layer k and k + 1 at nth time step, which can be calculated by F n k+1,I = −κ k,I T n k+1,I − T n k,I 1 where κ k,I is the thermal conductivity of the subsurface layer k of the patch I. The surface energy balance fluxes physically determine the conductive heat flux at the top of each patch as follows: where K is the number of subsurface layers. The bottom boundary condition of each patch uses prescribed interior temperatures for roofs and walls and zero temperature gradient for roads following [64].

Field Measurements
Field measurements conducted in a real urban environment are used to validate the performance of the MUSE model in simulating spatial heterogeneities of surface temperatures and surface energy balance fluxes. Field measurements were conducted at the Kongju National University campus in South Korea on a hot summer day of 19 August 2016 (Figure 5a). The study area consists of several buildings of apartments, campus buildings, and low residential houses. A micrometeorological tower has been installed at the rooftop of a campus building (36.4717 • N, 127.1426 • E), from which various meteorological and radiative variables have been obtained such as winds, air temperature, specific humidity, shortwave, and longwave radiative fluxes, and turbulent sensible and latent fluxes [69]. The surface temperatures were measured at six points on the building surfaces, with two horizontal surfaces (a roof and a road) and four vertical surfaces at different directions (east, west, south, and north) using a dual-point laser infrared thermometer (SK-8140, Sato Keiryoki Mfg. Co., Ltd.). The infrared thermometer has a measurement accuracy of ±2 • C. The roof of the building is made of concrete covered with green waterproof paint, and the wall surfaces are covered with red bricks. The road surface is made of asphalt.

Configuration of the Model
The simulation domain of the MUSE model configures with a 138 × 112 × 18 grid, which covers 414 × 336 m 2 horizontally and 54 m vertically with the same grid resolution of 3 m (Figure 5b). The resultant total number of active patches is 24,323 in the simulation domain. The physical property parameters of each patch have been assigned according to the surface features of the real urban area. The aerodynamic roughness lengths for momentum are 0.05 m for roof and road patches, and the roughness lengths for heat are 0.0005 m for roof patches and 0.005 m for road patches. The radiative and thermal parameters of roof and wall patches are assigned for concrete and red brick, and those of road patches are for asphalt [64,70]. Table 3 summarizes the physical property parameters assigned for the real urban simulation. The roof, wall, and road patches are 0.5, 0.4, and 1.0 m in depth, configuring with ten sublayers at each patch. Each surface's depth was determined based on the previous simulation studies that showed reasonable comparison against measurements (e.g., [64]). The tower measurements are used to provide accurate meteorological forcing for validation of the MUSE model. The global shortwave and longwave radiative fluxes (S dir↓ , S di f ↓ , and L di f ↓ ) from a four-component radiometer (CNR4, Campbell Scientific Inc., Cache Valley, UT, USA) deployed at the meteorological tower are applied homogeneously to the active patches, in which the measured global shortwave radiative flux is partitioned to the direct and the diffuse radiative fluxes following the empirical method [43,71]. The wind speed (U a ) and air temperature (T a ) adjacent to the active patches are obtained with interpolation from the wind speed and air temperature profiles computed using the flux measurements from the three-dimensional sonic anemometer (CSAT3, Campbell Scientific Inc., USA) and the Monin-Obukhov similarity relations as follows: where U(z) and T(z) are the wind speed and air temperature at height z, respectively, T obs is the air temperature measured at the height z obs (=33 m a.g.l.). u * , T * , and L are the friction velocity, the surface-layer temperature scale, and the Obukhov length, which were measured directly from the flux tower. The aerodynamic roughness length is set to 1.06 m following the simple formulation of Grimmond and Oke [72] from the analysis of the site measurement [69], and the integrated stability functions for momentum (ψ m ) and heat (ψ h ) are given in Högström's work ( [73] and [74], respectively).
The estimated wind speed and air temperature profiles are then horizontally interpolated to the model grid homogeneously as a meteorological forcing. In this study, the CFD model developed by Kim et al. [75] is used to force the MUSE model. The model is simultaneously configured with a 138 × 112 × 19 grid using the same grid resolution of the MUSE model based on the Arakawa-C type grid system, which is adopted frequently in many CFD and LES models (e.g., [3,65]). These grid configurations make the MUSE model directly compatible with the micrometeorological models and available "online" and "offline" mode computation. Here, the CFD model ran with the estimated wind speed profile under the assumption of neutral atmospheric stratification, and the resultant wind fields were used as meteorological forcing of the MUSE model. The wind fields simulated by the CFD model are more realistic in spatial distribution over the simulation domain than the homogeneously distributed wind speed profile.
The MUSE validation simulation was conducted for 28 h from 00 LST 19 August to 04 LST 20 August 2016. The surface temperatures of every active patch were initialized from the spin-up simulation of the previous 18 days. The validation simulation marched with the time step of 5 sec and took approximately 25 min using a single CPU. Figure 6 shows diurnal variations of the meteorological forcing applied for the MUSE validation simulation. The air temperatures measured at 33 m a.g.l. ranged diurnally from 24.6 • C to 33.9 • C, and the nocturnal temperatures stayed above 26.2 • C (a tropical night) on 20 August 2016. The wind speed measured at 33 m a.g.l. remained quite low throughout the day, mostly ranging 0.5-1.0 m s −1 (Figure 6a). The low wind speed was found frequently in forming the surface layer on a hot summer day. The downward shortwave radiation reached approximately 825 W m −2 around noon, but it was influenced by clouds observed during the measurement on the day (Figure 6b). The wind speed and air temperature profiles estimated from the measurements show characteristic diurnal variations within the surface layer on the day (Figure 6c,d). The wind speeds are quite low near the ground throughout the day, which gradually increases with height. The air temperature profiles are relatively homogeneous in the vertical direction during the night, which indicates that the atmosphere remained near neutral or weakly unstable.     (Figure 8a), while the nonsymmetrical street canyon has the two asymmetric building heights of 10 and 150 m, the canyon width of 30 m, and the canyon length of 45 m (Figure 8b). The MUSE model computes the patch view factors with a grid spacing of 1 m for both the street canyons. Meanwhile, the formulation of Johnson and Watson [59] gives the analytical sky view factors for validation. Figure 8 compares the numerical and analytical patch-sky view factors for the different geometric canyons. The sky view factors on the street canyon gradually decrease across the street canyon and increase along the street canyon with the distance from the canyon center. The patch-sky view factors at the higher building side are higher for the nonsymmetrical street canyon (Figure 8b). The MUSE model agrees well with the analytical patch-sky view factors at all the patches on both the symmetrical and nonsymmetrical street canyons. The formulations of Lee et al. [63] (Equations (12) and (15) in their paper) showed excellent agreement with the relative errors of less than 0.2% in the average sky view factor for the ideal street canyons, but they are applicable only for computing the sky view factors. The MUSE model shows a better agreement with relative errors of only less than 0.1% ( Table 4). The MUSE model's high accuracy stems from the formulations derived based on analytical relations [63]. Note that because the method's error can increase with patch size and short distance, the computation of the view factors for two adjacent patches was calculated by dividing the adjacent patches into 3 × 3 subpatches.  (12) in Lee et al. [63] is also compared.

East-West Direction North-South Direction Whole Area
Symmetric canyon    Figure 9 presents the spatial distribution of the patch-sky view factors computed for the real urban area. The patch-sky view factors at horizontal patches have values ranging from 0 to 1, whereas the patch-sky view factors at vertical wall patches have values ranging from 0 up to 0.5. It shows that the patch-sky view factors vary significantly over the domain by representing the complicated building morphology. The patch-sky view factor is low in dense building areas while it is relatively high on the roof patches and open ground areas. Tall buildings influence the surrounding low-rise buildings by reducing the sky view factors. As the MUSE model can calculate the patch view factors on an analytical level of accuracy for the idealized urban geometries (Table 4), the computation uncertainty for the real urban area is only associated with the representation accuracy of buildings' geometries, showing the capability of the MUSE model in determining the patch view factors for real urban environments. The effective albedo, which is defined by the ratio of outgoing shortwave radiation to incident shortwave radiation, is an appropriate parameter for validation of the shortwave radiative processes within urban canyons, including radiative trapping effects discussed in Section 2.2.3. The experimental study of Aida [76] is used for the validation, which has been applied in many previous studies (e.g., [40,67,77,78]). In the experimental study, ideal canyon structures ("model") of the northsouth canyon model ("model 1") and the east-west canyon model ("model 2"), and the cubed array canyon model ("model 3") were set up using cubed concrete blocks with height and width of 0.15 m. The effective albedo was measured at 0.45 m high using pyranometers. The MUSE model is applied for the experimental data obtained on a summer day (15 June 1978) and a winter day (3 December 1977) to the three different canyon geometries, which are configured with the same size of blocks as in Aida [76] experiments. The block surface albedos of wall and ground surfaces are assigned by 0.4 following Aida [76], but the rooftop and road surface albedos are assigned to vary diurnally according to an empirical correction relation of the albedo following Sievers and Zdunkowski [77]. It was found that the surface albedo varied with the altitude of the sun due to the angular dependency of surface reflectivity in the experiment for the flat surface model ("model 0") (Figure 10a,b). However, the MUSE model can not capture the angular dependency because it assumes Lambertian isotropic surface reflection. Therefore, the empirical corrections of the flat surface albedos were applied for validation purposes. The downward shortwave radiation and its partitioning into the direct and diffuse shortwave radiation are calculated based on empirical formulas in Overby et al. [43] with clear-day atmospheric turbidity conditions in Monteith and Unsworth [71].
The temporal variations of the effective albedos simulated for the three canyon geometries The effective albedo, which is defined by the ratio of outgoing shortwave radiation to incident shortwave radiation, is an appropriate parameter for validation of the shortwave radiative processes within urban canyons, including radiative trapping effects discussed in Section 2.2.3. The experimental study of Aida [76] is used for the validation, which has been applied in many previous studies (e.g., [40,67,77,78]). In the experimental study, ideal canyon structures ("model") of the north-south canyon model ("model 1") and the east-west canyon model ("model 2"), and the cubed array canyon model ("model 3") were set up using cubed concrete blocks with height and width of 0.15 m. The effective albedo was measured at 0.45 m high using pyranometers. The MUSE model is applied for the experimental data obtained on a summer day (15 June 1978) and a winter day (3 December 1977) to the three different canyon geometries, which are configured with the same size of blocks as in Aida [76] experiments. The block surface albedos of wall and ground surfaces are assigned by 0.4 following Aida [76], but the rooftop and road surface albedos are assigned to vary diurnally according to an empirical correction relation of the albedo following Sievers and Zdunkowski [77]. It was found that the surface albedo varied with the altitude of the sun due to the angular dependency of surface reflectivity in the experiment for the flat surface model ("model 0") (Figure 10a,b). However, the MUSE model can not capture the angular dependency because it assumes Lambertian isotropic surface reflection. Therefore, the empirical corrections of the flat surface albedos were applied for validation purposes.
The downward shortwave radiation and its partitioning into the direct and diffuse shortwave radiation are calculated based on empirical formulas in Overby et al. [43] with clear-day atmospheric turbidity conditions in Monteith and Unsworth [71]. The temporal variations of the effective albedos simulated for the three canyon geometries compare well with the measurements, showing that the MUSE model captures the characteristic temporal variations measured at the three canyon geometries in different seasons (Figure 10c-h). It has the temporal correlation coefficients of 0.69-0.99 for the summer day and 0.85-0.92 for the winter day. The relatively low correlation found in model 3 may attribute to the fact that the variation of the measured effective albedos is flat during the daytime. Meanwhile, the mean absolute errors of the simulated effective albedos are 0.002-0.016 for the summer day and 0.009-0.014 for the winter day, which corresponds to only 0.6-5.5% errors against the measurements (Table 5). It implies that the model errors are almost within the measurement accuracy (~5%).

The Real Urban Simulation
The model's performance in simulating surface temperatures and surface energy balance fluxes in real urban environments was tested against the measurements obtained from the Kongju National University area on a hot summer day. Figure 11 demonstrates spatial distributions of the simulated surface temperatures at different times of the day. It shows that the MUSE model can capture the evolution of spatial and temporal variations of the surface temperatures on each patch level. The simulated surface temperatures range spatially by 34 ± 4 • C at 09 LST, 49 ± 12 • C at 12 LST, 49 ± 10 • C at 15 LST, and 39 ± 5 • C at 18 LST. The surface temperatures and their spatial heterogeneities are relatively high around noon, as was discussed in Kim et al. [75]. Relatively lower surface temperatures are found at self-shaded walls and densely built-up areas, whereas the open patches such as rooftops and broad roads have relatively higher surface temperatures. These temporal and spatial variations are closely correlated with the solar radiation amount reaching the patches, suggesting the importance of accurate representation of the radiative processes in the model. Figure 12 compares the diurnal variations of the simulated surface temperatures against the field measurements at the roof, road, and four walls surfaces. The wall surface temperature has a characteristic diurnal variation according to its facing orientation (Figure 12a-d). The sun begins to heat the east wall first in the morning, and it heats the south and west walls in order as time goes on. The north wall has the smallest diurnal variation as expected. The roof and road surface temperatures have relatively higher surface temperatures due to their openness (Figure 12e,f). The model tends to slightly underestimate the wall temperatures while it overestimates the peak temperatures at the roof and road surfaces. The model measurement discrepancies might attribute to many aspects associated with the physical processes such as radiation, turbulent heat exchange, and subsurface heat storage, but they are not quantitatively discernible in this study. The surface temperatures simulated by the MUSE model compare well with the measurements, capturing the characteristic temporal patterns at the six different patch surfaces with correlation coefficients in the range of 0.94-0.98. The statistical evaluation results show that the mean bias errors (MBE) are only less than 3 • C and the root mean square errors (RMSE) are in the range of approximately 2-4 • C (Table 6). simulated surface temperatures range spatially by 34 ± 4 °C at 09 LST, 49 ± 12 °C at 12 LST, 49 ± 10 °C at 15 LST, and 39 ± 5 °C at 18 LST. The surface temperatures and their spatial heterogeneities are relatively high around noon, as was discussed in Kim et al. [75]. Relatively lower surface temperatures are found at self-shaded walls and densely built-up areas, whereas the open patches such as rooftops and broad roads have relatively higher surface temperatures. These temporal and spatial variations are closely correlated with the solar radiation amount reaching the patches, suggesting the importance of accurate representation of the radiative processes in the model.  Figure 12 compares the diurnal variations of the simulated surface temperatures against the field measurements at the roof, road, and four walls surfaces. The wall surface temperature has a characteristic diurnal variation according to its facing orientation (Figure 12a-d). The sun begins to heat the east wall first in the morning, and it heats the south and west walls in order as time goes on. The north wall has the smallest diurnal variation as expected. The roof and road surface temperatures have relatively higher surface temperatures due to their openness (Figure 12e,f). The model tends to slightly underestimate the wall temperatures while it overestimates the peak temperatures at the roof and road surfaces. The model measurement discrepancies might attribute to many aspects associated with the physical processes such as radiation, turbulent heat exchange, and subsurface heat storage,  (Table 6).     Figure 13 compares the diurnal variations of the simulated surface energy balance fluxes against the flux tower measurements at the urban site, and the statistical evaluation results are given in The statistical analysis shows that the model biases are systematic from the ratio of RMSE s to RMSE by 85% and 61% both in net radiation and storage heat flux, respectively. The temporal correlation coefficients are high-0.99 and 0.96 for net radiation and storage heat flux, respectively. These results imply that the model measurement discrepancies might be associated with static surface physical properties over the domain, such as surface albedo and thermal heat capacity. Despite the discrepancies, the MUSE model reproduces the surface energy balance over the area. The measured and simulated ratios of the storage heat flux to net radiation (G/R*) are about 0.72 and 0.69 during the daytime (0800-1600 LST) and 1.06 and 1.13 at night (2000-0200 LST), respectively. The relatively high G/R* is a characteristic feature of the heatwave period associated with weak winds in the urban area. Overall the MUSE model simulates diurnal variations of the surface fluxes and the relative energy balance measured in the urban area.

Summary and Conclusions
The MUSE model that is capable of representing the urban physical processes and their interactions in complex urban surfaces, and that can be compatible easily with general Cartesian grid microscale atmospheric models, was developed, which was coded in Fortran. The model represents complex urban surface morphology using three-dimensional patch surfaces, and it solves the urban physical processes of shortwave and longwave radiative transfers, turbulent exchanges of momentum and heat, and conductive heat transfer into urban subsurfaces. A shadow submodel is embedded to calculate time-varying shadow distributions in general urban environments, and the analytically-based numerical formulation is applied to estimate three-dimensional patch view factor distributions. The Monin-Obukhov similarity relation and the empirical formulation is used for turbulent heat exchanges on the horizontal (roofs and roads) and vertical wall surfaces, respectively. The conductive heat transfer at every patch is calculated based on Fourier's law. The MUSE model solves all of the urban physical processes interactively during the time integration and applies computationally efficient treatments for real urban applications.
We evaluated the MUSE model in terms of all of the model's physical processes and the full capability for use in a real urban environment. The embedded physical processes' validation results showed that the model could compute three-dimensional building shadows, patch view factors, and the urban radiative transfers for ideal urban geometries on a high accuracy level comparable to analytic solutions or measurement errors. The model's full capability was validated against the field measurements in simulating three-dimensional surface temperatures and surface energy balance fluxes at the Kongju National University campus area on a 3-m grid resolution. It first demonstrated that the model could reproduce the three-dimensional distributions of building shadows and patch view factors reflecting the urban area's morphological features. The simulation results showed that the surface thermal states changed significantly in space and time, and the model captured the thermal heterogeneities over the campus area. The MUSE model's realistic reproduction might be attributed to the interactive parameterization of the urban physical processes in real urban environments. Overall, the validation results demonstrate potential capabilities of the MUSE model for use in real urban simulations, implying that the assumption of homogeneous surface forcing found frequently in microscale atmospheric modeling should be reconsidered for realistic urban flow and dispersion simulations.
Because the model grid structure was designed to be compatible directly with general Cartesian grid microscale atmospheric models, the MUSE model can be easily used with CFDs and LESs models. The validation simulation in this study showed the capability in providing the mechanical forcing of every patch over the domain. The integrated modeling frame will be useful in interpreting various urban atmospheric and environmental problems (e.g., ventilation, heatwave, air quality) in real urban environments (e.g., [75,79]). For example, Kim et al. [75] examined the ventilation and thermal discomfort in a highly urbanized area using the CFD and MUSE coupled model.
As most CFD and LES models that have been used in urban microscale modeling have been developed based on a Cartesian coordinate, all the buildings in MUSE are thus represented by E-W and N-S parallel patches accordingly in this study, and the current realization of the MUSE formulation is limited to only horizontal and vertical surfaces. Realistic representation of slope surfaces and associated physical processes will be necessary for more realistic simulation, and the MUSE formulation given in this study can be applied to slope surfaces as well in calculation of shadow, view factor, and associated physical processes. The urban vegetation and associated hydrological processes are also important meteorological and environmental ingredients (e.g., [64,79,80]). Further development of the MUSE model for a better representation of real urban surfaces is currently underway.