Simulation of Spectral Albedo and Bidirectional Reflectance over Snow-Covered Urban Canyon: Model Development and Factor Analysis

Abstract

A critical comprehension of the impact of snow cover on urban bidirectional reflectance is pivotal for precise assessments of energy budgets, radiative forcing, and urban climate change. This study develops a numerical model that employs the Monte Carlo ray-tracing technique and a snow anisotropic reflectance model (ART) to simulate spectral albedo and bidirectional reflectance, accounting for urban structure and snow anisotropy. Validation using three flat surfaces and MODIS data (snow-free, fresh snow, and melting snow scenarios) revealed minimal errors: the maximum domain-averaged BRDF bias was 0.01% for flat surfaces, and the overall model-MODIS deviation was less than 0.05. The model’s performance confirmed its accuracy in reproducing the reflectance spectrum. A thorough investigation of key factors affecting bidirectional reflectance in snow-covered urban canyons ensued, with snow coverage found to be the dominant influence. Urban coverage, building height, and soot pollutant concentration significantly impact visible and infrared reflectance, while snow grain size has the greatest effect on shortwave infrared. The bidirectional reflectance at backward scattering angles (0.5–0.6) at 645 nm is lower than forward scattering (around 0.8) in the principal plane as snow grain size increases. These findings contribute to a deeper understanding of snow-covered urban canyons’ reflectance characteristics and facilitate the quantification of radiation interactions, cloud-snow discrimination, and satellite-based retrieval of aerosol and snow parameters.

1. Introduction

The world has undergone a dramatic process of urbanization, with half of the world’s population living in urban areas [1] and more than 600 million people now living in winter cities located above 45° worldwide [2]. Man-made urban surfaces (e.g., concrete pavements, buildings, and asphalt roads) generally have higher surface reflectance than natural surfaces (e.g., grass, forest, and wetlands). Long-term snow cover over urban canyons during wintertime further increases urban surface reflectance and redistributes the radiation budget [3,4], and urbanization-induced emissions of anthropogenic pollutants and enhanced anthropogenic heating together lead to varying degrees of impact on local weather, environment, and climate [5]. Therefore, it is important to clarify the influence of snow cover on the reflectivity of urban canyons for a better understanding of the urban climate system.

In addition, snow-covered surfaces are also a great challenge for aerosol remote sensing [6]. Snow cover over urban canyons creates strong heterogeneity and different reflectivity characteristics compared to snow-free urban canyons and flat snow surfaces, and processes such as pollutant deposition, snow grain aging, and snow melting further change their directional reflectivity characteristics [3,7]. Separating and explicitly describing the aerosol contributions from the top-of-atmosphere (TOA) radiance measured by satellite-based instruments is at the core of aerosol remote sensing, but it is still challenging because the contribution of highly reflective snow-covered surfaces is much higher than that of aerosol, which could easily lead to unacceptable inversion bias [8]. Therefore, a thorough understanding of the influence of snow cover on the reflectance of urban canyons is essential for further improvements in aerosol remote sensing in winter cities.

The bidirectional reflectance distribution function (BRDF) is a critical factor in aerosol remote sensing, radiative forcing, and albedo feedback [9,10]. The BRDF model has the potential to reduce errors in mapping snow-covered areas, separating snow from clouds, and retrieving aerosol optical depth (AOD) [11]. For example, Ji et al. [12] collected and compared the differences in BRDF between pure snow, pure ice, and ice-snow to support snow/ice/ice-snow classification and monitoring. Gatebe and King [13] established an airborne BRDF database of snow, cloud, and other major surface types over a 30-year period (1984–2014) to support snow/cloud discrimination and various remote sensing applications. Tian et al. [14] and Lin et al. [15] reported the potential use of BRDF in improving the quality of AOD retrieval over urban areas.

Efforts in remote sensing of aerosols over the cryosphere using snow BRDF have been reported in many previous studies [6,16]. For example, Shi et al. [17] used the dual-angle observation of AATSR to perform AOD inversion in the Arctic region, and the correlation coefficients between the inversion results and the ground-based observations were 0.75 and 0.70, respectively. Mei et al. [8] improved the XBAER inversion algorithm to achieve AOD inversion over the cryosphere, which assumes that the surface reflectance can be characterized by a one-dimensional linear equation, and the regression coefficients are stored in a pre-calculated database based on a large number of observations. Swain et al. [18] retrieved the AOD over snow-covered regions of the Arctic using passive remote sensing from space, and the surface reflectance was parameterized by a two-parameter snow BRDF model [19].

Research on snow BRDF models to support AOD inversion in polar regions has made good progress, but little has been reported on AOD inversion over snow-covered urban canyons. Establishing an accurate BRDF database to support AOD inversion, as Deep Blue [20], Dark Target [21], and XBAER algorithms [8] have previously done, is a challenging task. Compared with polar and alpine mountain environments, snow in winter cities is more often affected by aerosol deposition contamination, as well as grain aging and melting, resulting in continuous changes in reflectance properties [22]. Zhang et al. [23] reported that the deposition of black carbon and sandy aerosol particles led to a 20% decrease in snow albedo in the southeastern region of the Qinghai-Tibetan Plateau. Wang et al. [24] found that snow albedo decreased by 2–6% due to BC contamination with 100–500 ng/g⁻¹ BC in the snowpack. Thus, the surface reflectance database constructed under the assumption of short-term invariance cannot accurately characterize the time series variations of the snow-covered urban canyon BRDF, resulting in the lack of winter AOD products over high latitude cities. The lack of in-depth understanding of the effect of snow on the spectral albedo and BRDF characteristics over the urban canyon is considered to be the primary factor.

Motivated by the need for BRDF characteristics over snow-covered urban canyons, we attempt to develop a new snow-covered urban bidirectional reflectance model to comprehensively evaluate the effect of snow cover on urban reflectance. This model is established based on the Monte Carlo ray tracing method, which has been widely used and validated to simulate directional reflectance over rugged terrain [25,26,27] and terrestrial snow [28,29,30]. The asymptotic radiative transfer (ART) model is applied to simulate the bidirectional reflectivity characteristics of snow, since the ART model has been extensively validated to retrieve the snow grain size, pollutant concentration, and snow albedo [31,32,33]. This study is organized as follows. Section 2 gives detailed descriptions of the data and methods used, including Monte Carlo ray tracing and the ART model. Section 3 presents the model validation, and Section 4 gives further assessments of the influence of key factors on urban BRDF. Section 5 gives a brief conclusion and prospects.

2. Materials and Methods

2.1. Site Description

Harbin, often referred to as the “Ice City” of China, is the provincial capital of Heilongjiang and is situated in the northernmost part of the Northeast China Plain (Figure 1a). With a temperate continental monsoon climate featuring approximately four months of snow cover [34,35], Harbin is a megacity spanning over 10,200 square kilometers and home to approximately 10.7 million permanent residents [36]. Given its prolonged snow cover and moderate air pollution [37,38], we selected this region (125°42′E–130°10′E, 44°04′N–46°40′N) to validate our numerical model. Figure 1c–e shows the corresponding spectral albedo map on 13 October 2019 (autumn season). The domain-averaged albedos at 645 nm, 469 nm, and 1640 nm are 0.098, 0.053, and 0.220, respectively. In Figure 1e, the Songhua River and urban areas stand out with lower reflectance in the shortwave infrared bands compared to the surrounding vegetation.

2.2. Geometric Modeling of Urban Canyons

Considering the variations in urban architecture, such as commercial centers, residential zones, and factories, we design an urban canyon model with a 25% building coverage to accurately represent the urban layout and ensure numerical simulations are representative. The urban canyon, measuring 2.5 km by 2.5 km, is divided into a 50 × 50 m grid in both the x and y directions. Building heights are assumed to follow a normal distribution, with an average of 50 m and a standard deviation of 0.6. The three-dimensional representation of the urban canyon is displayed in Figure 2.

The geometric modeling defaults and their variation ranges for urban canyons and snow are detailed in

Table 1. Input parameters of the numeric model.

Parameter	Unit	Range	Default
building coverage	-	10–50%	25%
building height	meter	10–90	50
length of x	km	-	2.5
length of y	km	-	2.5
snow coverage	-	20–100%	80%
snow grain size	µm	1000–5000	3600
pollutant concentration	ppb	${10}^2$–${10}^6$	500
zenith angle	degree	40–80	70
azimuth angle	degree	0–180	0

. The default parameters include a building coverage (BC) of 25%, a building height (BH) of 50 m, a solar zenith angle (SZA) of 70°, and a relative azimuth angle of 0°. To investigate the impact of varying building structure and snow conditions on spectral albedo and BRDF properties, we consider the following ranges: building coverage from 10% to 50%, building height between 10 m and 90 m with a standard deviation of 0.6, solar zenith angle ranging from 40° to 80° in increments of 10°, and relative azimuth angle varying from 0° to 180°.

For snow, we focus on snow coverage (SC), snow grain size (SGS), and soot pollutant concentration (SPC). The initial settings are 80% snow coverage, 3600 µm snow grain size, and a soot pollutant concentration of 500 ppb. In our analysis, snow coverage varies between 20% and 100%, snow grain size ranges from 100 µm to 5000 µm, and soot pollutant concentrations span from 10² ppb to 10⁶ ppb.

2.3. Satellite Albedo Product

The Terra MOD09A1 and Aqua MYD09A1 Version 6 Moderate Resolution Imaging Spectroradiometer (MODIS) products provide surface spectral albedo estimates for Bands 1 through 7 (as listed in

Table 2. Band information of the MODIS albedo product.

Band	Wavelength	Spectra
Band 1	620–670 nm	Red
Band 2	841–876 nm	NIR
Band 3	459–479 nm	Blue
Band 4	545–565 nm	Green
Band 5	1230–1250 nm	SWIR
Band 6	1628–1652 nm	SWIR
Band 7	2105–2155 nm	SWIR

). These data have been meticulously corrected for atmospheric influences, such as gases, aerosols, and Rayleigh scattering. For each pixel, the most suitable observation within an eight-day composite period is chosen based on cloud cover and solar zenith angle. In cases where multiple acquisitions meet the criteria, the pixel with the lowest Band 3 (blue) value is selected. This choice represents the optimal Level 2G (L2G) observation, which maximizes coverage, avoids clouds and cloud shadows, and minimizes aerosol effects. The Collection 6 surface reflectance product exhibits an accuracy of at least 0.005 + 0.05 SR across all bands, where SR denotes surface reflectance [39].

The MODIS spectral albedo is derived from an operational algorithm that integrates the RossThick–LiSparseReciprocal (RTLSR) semi-empirical model for inversion [40,41]. This model decomposes the bidirectional reflectance distribution function (BRDF) into a weighted sum of an isotropic component ($f_{iso}$) and two distinct kernels. One, $K_{vol}$, is based on volume scattering radiative transfer models, while the other, $K_{geo}$, is informed by surface scattering and geometric shadow casting theory [42]. The operational product computes BRDF parameters ($f_{iso}$, $f_{vol}$, $f_{geo}$) that are spectrally dependent weights assigned to each kernel for calculating the overall anisotropic reflectance. Figure 3 illustrates the viewing geometry and planes, including the principal plane (PP), cross plane (CP), and almucantar plane (ALM).

$$\begin{array}{cc}\hfill R({\theta }_s,{\theta }_v,\phi ,\lambda )=& f_{iso}\left(\lambda \right)\hfill \\ \hfill & +f_{geo}\left(\lambda \right)K_{geo}({\theta }_s,{\theta }_v,\phi )\hfill \\ \hfill & +f_{vol}\left(\lambda \right)K_{vol}({\theta }_s,{\theta }_v,\phi ),\hfill \end{array}$$

(1)

where ${\theta }_s$, ${\theta }_v$, $\phi$, and $\lambda$ represent the solar zenith angle, view zenith angle, relative azimuth angle, and wavelength, respectively. $K_{geo}({\theta }_s,{\theta }_v,\phi )$ and $K_{vol}({\theta }_s,{\theta }_v,\phi )$ are the geometric optical kernel and the volume scattering kernel, respectively. $f_{iso}$, $f_{geo}$, and $f_{vol}$ represent the weight coefficients of isotropical kernel (equal to 1), geometric optical kernel, and scattering kernel with wavelength $\lambda$, respectively. The RossThick volume scattering kernel and the LiSparseReciprocal geometric optical kernel are formulated as follows:

$$\begin{array}{cc}\hfill K_{geo}({\theta }_s,{\theta }_v,\phi ) =& \mathrm{O}({\theta }_{\mathrm{s}},{\theta }_v,\phi )-sec{\theta }_{\mathrm{s}}-sec{\theta }_v+\hfill \\ \hfill & \frac{1}{2}(1+cos\xi )sec{\theta }_{v}sec{\theta }_{\mathrm{s}}\hfill \end{array}$$

(2)

$$\begin{array}{c}\hfill K_{vol}({\theta }_s,{\theta }_v,\phi )=\frac{(\pi /2-\xi )cos\xi +sin\xi }{cos{\theta }_s+cos{\theta }_v}-\frac{\pi }{4},\end{array}$$

(3)

where $\mathrm{O}({\theta }_{\mathrm{s}},{\theta }_v,\phi )$ is the overlap function of the viewing and sunlit shadows, and $\xi$ refers to the phase angle regarding the viewing and illuminating geometries [42].

To facilitate urban pixel extraction, a 500-meter land cover dataset is utilized. The MCD12Q1 Version 6 data product, which offers annual land cover classifications, is derived from MODIS Terra and Aqua reflectance data through supervised classification methods [43]. The International Geosphere-Biosphere Programme (IGBP) legend is employed to identify urban areas. This classification scheme consists of 18 categories, with urban and built-up lands designated as Class 13. This classification is based on the presence of at least 30% impervious surface area, including buildings, asphalt, and vehicles [44,45].

2.4. Snow ART Model

Snow is primarily composed of ice crystals. Upon interaction with sunlight, these crystals predominantly reflect photons due to their icy surfaces, which then undergo reflection and refraction. The ART model, based on the concept of fractal particles, estimates snow reflectance by considering the influence of soot pollution and the size of snow grains. The bidirectional reflectance of snow is approximated using an exponential function [19,46]:

$$r_{\mathrm{so}}=r_{0}exp\left(\frac{-\alpha K_{\mathrm{s}}{K}_{\mathrm{v}}}{r_0}\right)$$

(4)

$$K_{\mathrm{s}}=\frac{3}{7}(1+2cos{\theta }_{\mathrm{s}}^{\prime })$$

(5)

$$K_{\mathrm{v}}=\frac{3}{7}(1+2cos{\theta }_{\mathrm{v}}^{\prime })$$

(6)

$$r_0=\frac{a+b(cos{\theta \prime }_{\mathrm{s}}+cos{\theta \prime }_{\mathrm{v}})+c·cos{\theta \prime }_{\mathrm{s}}cos{\theta \prime }_{\mathrm{v}}+p}{4(cos{\theta \prime }_{\mathrm{s}}+cos{\theta \prime }_{\mathrm{v}})}$$

(7)

$$p=11.1exp(-0.087\Theta )+1.1exp(-0.014\Theta )$$

(8)

$$cos\Theta =-cos{\theta }_{\mathrm{s}}^{\prime }cos{\theta }_{\mathrm{v}}^{\prime }+sin{\theta }_{\mathrm{s}}^{\prime }sin{\theta }_{\mathrm{v}}^{\prime }cos\left({\phi }_s^{\prime }-{\phi }_{\nu }^{\prime }\right)$$

(9)

$$\alpha =A\sqrt{\frac{4\pi a_{\mathrm{ef}}}{\lambda }(k+C_{\mathrm{soot}})},$$

(10)

where ${\theta \prime }_{\mathrm{s}}$ and ${\theta \prime }_{\mathrm{v}}$ are relative solar zenith angle and relative view zenith angle, respectively. $K_{\mathrm{s}}$ and $K_{\mathrm{v}}$ are the escape functions for different relative solar zenith angles and view zenith angles, respectively. $r_0$ is the reflectance of zero absorption medium, while a, b, and c equal 1.247, 1.186, and 5.157, respectively. ${\phi }_s^{\prime }$ and ${\phi }_{\nu }^{\prime }$ are the relative solar and view azimuth angle, respectively. $\Theta$ is the scattering angle. A is equal to 5.8, while $a_{\mathrm{ef}}$ and $C_{\mathrm{soot}}$ are snow grain size and pollutant concentration, respectively. k is related to the ice refractive index [11,46]. The spectral reflectance of snow using default values of snow grain and soot pollutant concentration is shown in Figure 4.

2.5. Spectral Albedo Data for Urban Surfaces

The spectral reflectance data for urban surfaces were sourced from the United States Geological Survey (USGS) Spectral Library version 7 (https://www.usgs.gov/labs/spectroscopy-lab/science/spectral-library, accessed on 1 April 2024), which contains thousands of laboratory-measured spectra of minerals, vegetation, and man-made materials [47]. Covering wavelengths from 0.2 µm to 200 µm, the library facilitates the understanding of remotely sensed data by linking spectral features to chemical structures. In this research, we employed eight man-made materials (brushed aluminum, asphalt road, asphalt shingle, asphalt tar, concrete road, concrete surface, painted aluminum, and roofing felt) and three natural materials (snow melt, aspen leaves, and gray pine branches). Their respective spectral reflectances are displayed in Figure 5.

We assume that, under normal conditions, a building roof typically consists of brushed aluminum (30%), asphalt shingle (10%), asphalt tar (10%), painted aluminum (30%), and roofing felt (20%). The ground surface is composed of aspen leaves (10%), asphalt road (10%), concrete road (50%), and gray pine (30%). When snow covers both the roof and ground, the roof is primarily snow (80%), with brushed aluminum (10%) and painted aluminum (10%) contributing, while the ground surface consists of snow (70%), asphalt road (5%), gray pine (5%), and concrete road (20%). The building facade is characterized by a concrete surface.

2.6. Monte Carlo Ray-Tracing Method

The Monte Carlo ray-tracing method, a versatile approach for solving the radiative transfer equation in complex environments like urban canyons with diverse surface properties (roofs, ground, and facades) [48,49], accurately models the domain-averaged albedo and bidirectional reflectance distribution function (BRDF) over snow-covered canyons by simulating 3D photon transport. The core principle of this technique involves simulating photon interactions with surfaces, accounting for reflection. Note that the proposed model disregards absorption and scattering by atmospheric particles [50].

Beginning with a photon’s energy, position, and direction, the process entails the following steps: upon entering the domain from the solar incidence direction, the photon’s intersection with a surface is identified. Once the reflection point is known, the new photon direction is determined based on the local BRDF, and the photon’s weight is adjusted according to the surface albedo. This iterative process continues until the photon exits the physical boundaries or is fully absorbed. The overall albedo is calculated as the ratio of outgoing to incoming photon energy (as depicted in Figure 6), and the BRDF is derived accordingly:

$$\begin{array}{cc}\hfill \mathrm{BRDF}=& \frac{\pi M_r\left({\Omega }_r\right)}{M_0{\Omega }_{r}cos{\theta }_r}\hfill \\ \hfill =& \frac{\pi M_r\left({\Omega }_r\right)}{M_{0}cos{\theta }_{r}sin{\theta }_{r}d{\theta }_{r}d\phi }\hfill \end{array}$$

(11)

$$d\Omega =\frac{dA}{r^2}=sin\theta d\theta d\phi ,$$

(12)

where ${\Omega }_r$ represents the reflection direction, $M_r\left({\Omega }_r\right)$ is the energy of outgoing photons at ${\theta }_r$ and $\phi$, $M_0$ is the total incident energy, and $\Omega$ represents the solid angle. In the urban canyon, surfaces are primarily characterized by Lambertian properties, while snow cover exhibits anisotropic behavior according to the RTLSR model. For multi-component roof and ground surfaces, the surface albedo $\rho$ and bidirectional reflectance distribution function (BRDF) are calculated by considering the albedo weighting:

$$\begin{array}{c}\hfill \rho =\frac{\sum {\eta }_i{\rho }_i}{\sum {\rho }_i}\end{array}$$

(13)

$$\begin{array}{c}\hfill \mathrm{BRDF}=\frac{\sum {\eta }_i{\mathrm{BRDF}}_i}{\sum {\mathrm{BRDF}}_i},\end{array}$$

(14)

where ${\eta }_i$ represents the area share of the ith surface.

3. Model Validation

3.1. Validation over Flat Surfaces

To evaluate the model performance, we compare the Monte Carlo-generated albedo to the exact value for a flat Lambertian surface with an albedo of 0.2 under a solar zenith angle of 45°. The surface, measuring 2.5 km by 2.5 km, is divided into a 50 × 50 grid with 50-m intervals in both x and y directions. Figure 7 demonstrates the trend of domain-averaged relative error as a function of incident photon count. It is clearly shown that the simulation error decreases as the photon count increases. The modeled BRDF exhibits a relative error of 0.75% and 0.49% when using ${10}^7$ and ${10}^8$ photons, respectively. To strike a balance between accuracy and computational efficiency, we opt for ${10}^7$ photons for further simulations.

In Figure 8, a comparison of modeled BRDF distributions is shown with respect to established distributions for a flat Lambertian surface and two anisotropic surfaces. The analysis is based on MODIS Band 3 at a wavelength of 469 nm. The anisotropic surfaces are characterized by the RTLSR, model with $f_{iso}=0.091$, $f_{vol}=0.032$, and $f_{geo}=0.012$, and the ART model with a snow grain size of 500 µm and a soot pollutant concentration of 3600 ppb. The domain-averaged albedos for the Lambertian, RTLSR, and ART surfaces are 0.199, 0.078, and 0.914, respectively, with corresponding relative biases of −0.49%, −1.57%, and −0.13%. The maximum biases across all cases do not exceed 0.01%, while the maximum relative biases are within 2% for the Lambertian, 5% for the RTLSR, and 3% for the ART surfaces. These results indicate that the numerical model performs admirably in flat surface scenarios.

3.2. Validation with MODIS Albedo Observations

Figure 9 illustrates the time series of spectral albedo for Band 1 (red), Band 3 (blue), and Band 6 (SWIR). Winter exhibits significantly higher albedos compared to other seasons, persisting for over four months. The mean spectral albedo at Band 1 (645 nm) is approximately 0.4 during winter, with a minimal difference between the MODIS Terra and Aqua measurements. Band 3 (469 nm) shows more pronounced variations, while the albedo at Band 6 (1640 nm) exhibits a contrasting pattern, with low values in winter (around 0.1) and higher values in the rest of the year (approximately 0.2). These spectral patterns are closely linked to snow cover, as snow reflects strongly in visible and infrared bands but has a lower reflectance in the shortwave infrared band (as depicted in Figure 4 and Figure 5).

To thoroughly assess the model’s performance, three MODIS dataset scenarios were extracted: snow-free, fresh snow, and snow melt. These scenarios were employed to evaluate if the numerical model could accurately reproduce observed spectral albedos. The analysis focused on Harbin, and only pixels with IGBP land cover types of urban and built-up are considered. Figure 10 illustrates the spectral changes in MOD09, MYD09, and modeled albedo for each condition. On 11 November 2019, Harbin was in the early winter with no recent snowfall, showing a lower visible band albedo compared to infrared and shortwave bands. A snowfall on 28 December 2019 led to a significant increase in surface albedo across all bands. By 11 February 2020, after a week without fresh snow, the visible band albedo decreased to around 0.35–0.40, remaining notably higher than the snow-free period. The model accurately captures the spectral albedo variations across these scenarios. As seen in Figure 10d–f, the model’s deviation from observations is generally smaller than the discrepancies between the MOD09 and MYD09 data. The overall deviation between the model and observations does not exceed 0.05 for these scenarios, suggesting that the model effectively simulates spectral albedo and exhibits minimal deviation from real-world observations.

4. Results

This section presents a comprehensive analysis of spectral albedo and bidirectional reflectance in relation to a number of key factors, including building coverage (BC), building height (BH), snow coverage (SC), snow grain size (SGS), soot pollutant concentration (SPC), and solar zenith angle (SZA). The default input variables presented in Table 1 are employed in this analysis.

4.1. Building Coverage

The influence of building coverage on the snow-covered urban canyon’s albedo spectrum and bidirectional reflectance is showcased in Figure 11. Surface albedo across the visible and near-infrared wavelengths (400–1000 nm) varies from 0.50 to 0.70, as the building coverage increases from 10% to 50%. In contrast, the albedo values in the shortwave infrared band (1600–2200 nm) exhibit a more limited range of 0.10 to 0.15. As depicted in Figure 11a, the surface albedo over snow-covered urban canyons decreases as building coverage rises from 10% to 50%. This decline is mainly due to the increased number of reflections between roofs, ground, and faceted surfaces, which reduces the total energy carried by escaping photons. Figure 11b,c display the bidirectional reflectance distributions at MODIS Band 1 (645 nm) and Band 6 (1640 nm) in the hemispherical space. The average albedos for Band 1 and Band 6 are 0.578 and 0.131, respectively. Both polar plots exhibit angular variation, with higher reflectance at large backward scattering angles compared to the domain average, while forward scattering angles show lower values. In contrast to the strong backward scattering in fresh snow on flat surfaces (Figure 8i), urban structures significantly impact the directional reflectance characteristics.

The bidirectional reflectance changes with varying building coverages in ALM, PP, and CP are illustrated in Figure 11d–i. Reflectance variations at Band 1 are more pronounced than at Band 6. In the ALM, bidirectional reflectance is generally high at both forward and backward directions, with a slight dip between them. In the PP and CP, the impact of building coverage at a 0° view zenith angle (VZA) is minimal, and this effect increases with increasing VZA. The presence of buildings notably reduces backward scatter reflectance at higher VZA in the PP and CP. The PP’s forward scattering reflectance shows a decreasing trend, followed by an increase with VZA, while the CP consistently exhibits a decreasing trend.

4.2. Building Height

Figure 12 demonstrates the influence of building height on snow-covered urban canyon albedo and bidirectional reflectance. As the height rises from 10 m to 90 m, a gradual decrease in spectral albedo across MODIS Bands 1–7 is noticeable. This decline tapers off with escalating height. In comparison to Figure 11, similar albedo levels are observed across the MODIS bands. Surface albedos in the visible and infrared regions fluctuate between 0.50 and 0.70, while those in the shortwave infrared range from 0.08 to 0.15. The lower albedo at higher building heights can be primarily attributed to enhanced interactions between photons and various surfaces, such as roofs, ground, and facades. Figure 12b,c display the bidirectional reflectance at Bands 1 and 6 in the hemispheric space for a 90-meter tall building. The domain-averaged reflectance at Bands 1 and 6 are 0.532 and 0.121, respectively. The angular variation in bidirectional reflectance is more pronounced with a building coverage ratio of 25% and a height of 90 m compared to a coverage of 50% and a height of 50 m. This is evident from the bidirectional reflectance decrease at low VZA, while high VZA values exhibit more extreme values in both Band 1 and Band 6.

The bidirectional reflectance changes with varying building heights in the ALM, PP, and CP are shown in Figure 12d–i. As the building height increases, the ALM exhibits less variation in reflectance between 120° and 180° relative azimuth angles, while backward scattering directions show greater variability. The reduction in reflectance due to height increase is more significant at the 240° to 360° angle compared to 0° to 120°. In the PP, backward scattering directions have significant variability, while forward scattering shows less. This variation decreases with increasing height. The CP displays a gradual rise in reflectance from forward to backward scattering, a trend linked to specific urban canyon structures. Given that our simulations employ a fixed default urban model (except in Section 3.1), the default bidirectional reflectance distribution is primarily determined by urban canyon characteristics. Although variations in urban structures can alter the pattern, the focus on trends rather than absolute values or distributions enables a comprehensive understanding of the key factors affecting bidirectional reflectance over snow-covered urban canyons.

4.3. Snow Coverage

The impact of snow coverage on the albedo spectrum and bidirectional reflectance within snow-covered urban canyons is illustrated in Figure 13. As snow coverage escalates, the proportion of reflectance contributed by non-snow surfaces decreases accordingly. Snow coverage significantly affects spectral albedo, with increased coverage enhancing reflectance in the visible and near-infrared regions (400–1000 nm), while having a dampening effect on shortwave infrared bands (1600–2200 nm). The influence of varying snow coverage on albedo at around 1.45 µm is minimal, likely due to the low reflection properties of snow and other materials in the 1.4–1.5 µm range, as depicted in Figure 5.

From Figure 13b,c, it is evident that the bidirectional reflectance distribution pattern displays minimal variation. Figure 13d–i illustrate the bidirectional reflectance alterations with increasing snow coverage in the ALM, PP, and CP. While there are no significant impacts on the angular variation of bidirectional reflectance in ALM, PP, and CP at MODIS Band 3, Band 6 exhibits a different pattern. In the PP and CP, bidirectional reflectance at forward and backward scattering angles becomes less sensitive to snow coverage. At smaller SZAs below 30°, bidirectional reflectance is more susceptible to changes in snow coverage. As snow albedo is lower at Band 6 compared to materials like aluminum, roofing felt, asphalt roads, and concrete roads (as depicted in Figure 5), an increase in snow coverage leads to a decrease in reflectance. This overall rise in bidirectional reflectance can be attributed to the interaction between photons and surfaces with higher albedos.

4.4. Snow Grain Size

The spectral variation of surface albedo with respect to snow grain size is illustrated in Figure 14. It is evident that snow grain size exerts a relatively small impact on albedo in the visible and near-infrared regions (400–1000 nm), but its influence becomes pronounced in the shortwave infrared band (1600–2200 nm). The average bidirectional reflectance values for a snow grain size of 5000 µm at Band 1 and Band 6 are 0.554 and 0.112, respectively.

The angular variations of bidirectional reflectance in the ALM, PP, and CP are depicted in Figure 14d–i. In comparison to building coverage, height, and snow cover, snow grain size has a relatively smaller impact on the angular variation pattern at MODIS Band 3, particularly in ALM, PP, and CP. However, this influence is more pronounced at Band 6, and it decreases gradually with larger snow grain sizes. At Band 1, the backscattering reflectance for a 5000 µm input is slightly higher than that for 4000 µm in CP, but this distinction is not observed in Band 6. In Band 1, the backward scattering direction exhibits lower reflectance (0.5–0.6) than the forward scattering direction (around 0.8) in PP. Conversely, at Band 6, the difference in bidirectional reflectance between forward and backward scattering directions is notably reduced.

4.5. Soot Pollutant Concentration

In Figure 15a, the impact of pollutant concentration on the snow-covered urban canyon’s albedo spectrum is evident. The albedo, particularly in the visible and near-infrared (400–1000 nm) bands, is significantly affected, with an average albedo of approximately 0.6 for a concentration of 100 ppb, which doubles to around 645 nm when the concentration is ${10}^6$ ppb. However, this influence is negligible in the shortwave infrared (1600–2200 nm) range. The comparison in Figure 15b,c illustrates that, when a soot pollutant concentration of ${10}^6$ ppb is introduced, the bidirectional reflectance in Band 1 decreases notably across the entire hemisphere, as opposed to Figure 14b, where the difference in Band 6 remains minimal.

Figure 15d–i illustrate the bidirectional reflectance changes with varying pollutant concentrations in the ALM, PP, and CP. At Band 6, the reflectance remains virtually unchanged in all three regions as the concentration increases exponentially. In contrast, at Band 1, a noticeable decrease in bidirectional reflectance is observed when the concentration rises from ${10}^5$ ppb to ${10}^6$ ppb. The angular variation pattern at Band 1 shows no significant deviation with the exponential increase in pollutants. The increment in soot concentration leads to a diminishing difference in reflectance between the zenith and nearby VZAs.

4.6. Solar Zenith Angle

Figure 16 demonstrates the effect of SZA on albedo spectrum and bidirectional reflectance. As the SZA rises from 40° to 80°, a slight decrease in spectral reflectance is noted in the visible and infrared bands. The angular distribution of bidirectional reflectance across the hemisphere is shown in Figure 16b,c, with domain-averaged values of 0.598 for Band 1 and 0.141 for Band 6. The increase in SZA leads to enhanced forward scattering, primarily due to the strong forward scattering characteristic of snow at high incident angles [4,8,11].

The angular variation of bidirectional reflectance in ALM, PP, and CP is depicted in Figure 16d–i. The SZA’s impact on the ALM’s bidirectional reflectance distribution is relatively minor, with similar trends observed in both Band 1 and Band 6. In contrast, the PP exhibits significant changes in bidirectional reflectance with increasing SZA. At low angles (below 60°), backward scattering directions have notably lower reflectance compared to zenith and forward scattering. As the angle increases to 80°, forward and backward scattering directions experience a substantial boost in reflectance, while the reflectance at the zenith direction decreases. A similar pattern is observed in the CP, although the SWIR overall reflectance remains relatively constant. However, a significant rise in backward scattering reflectance is evident in the PP.

4.7. Influence of Intertwined Factors on Snow-Covered Urban Albedo

Two scenarios examining the combined effects of various factors on snow-covered urban albedo are explored. Figure 17a–c illustrate the spectral albedo changes with building height and snow coverage. Notably, higher snow coverage and lower building heights result in increased spectral albedo in the visible bands. In contrast, the shortwave infrared band exhibits significantly lower overall albedo, with a slight increase observed as building height and snow coverage decrease. Figure 17d–f highlight the impact of snow grain size and soot pollution concentration. Soot pollution obstructs upward reflection in snow-covered urban canyons. In the visible bands, albedo decreases with larger snow grain size and higher soot pollution. However, in the shortwave infrared band, the primary driver of albedo reduction is the increase in snow grain size, as the influence of soot concentration is negligible due to the low reflectance of snow in this band compared to the visible spectrum and other urban surfaces.

5. Discussions

Previous research has developed various bidirectional reflectance models for snow-covered environments like mountains and forests.

Table 3. A brief summary of bidirectional reflectance models for snow-covered regions.

No.	Surface	Method	Reference
1	sastrugi over rough snow	Photometric Roughness Model and a multiple scattering model	Leroux and Fily [51]
2	macroscopic rough snow	Monte Carlo ray-tracing method	Larue et al. [30]
3	small-scale rough snow	Ray-tracing method with photon recollision probability theory	Manninen et al. [53]
4	mountain	Monte Carlo ray-tracing method	Lee et al. [26], Liou et al. [52]
5	mountain	Monte Carlo ray-tracing method	Chu et al. [9]
6	mountain	RTLSR BRDF model	Bertoncini et al. [54]
7	forest	hybrid geometric optical-radiative transfer model	Chen et al. [11]
8	forest	photon recollision probability theory	Jääskeläinen and Manninen [55]
9	urban	Monte Carlo ray-tracing method	This Study

provides a concise summary of bidirectional reflectance models applicable to various snow-covered terrains, such as rough snow surfaces, mountainous regions, and forested areas. These models are founded on diverse methodologies, including Monte Carlo ray-tracing techniques [9,26,30,51,52,53], semi-empirical BRDF approaches [54], and geometric optics [11]. A consistent discovery across these studies is the significant variability in snow’s bidirectional reflectance, influenced by the complex interactions with rough surfaces, mountainous terrains, and forests. These interactions with the environment lead to alterations in the BRDF patterns to varying extents. However, a specific model for snow-covered urban canyons is lacking, hindering our understanding of snow’s impact on albedo, climate, and remote sensing applications such as aerosol, snow grain, and cloud detection in urban regions. To address this gap, we propose a new snow-covered urban canyon bidirectional reflectance model based on the Monte Carlo ray-tracing technique. The ART model characterizes snow reflectance, while other surfaces use the Lambert model. The findings will advance knowledge about snow’s effect on urban canyon anisotropy and inspire further research, including radiation interactions, urban reflectance sensitivity, cloud-snow discrimination, and satellite-based parameter retrieval (e.g., aerosol optical depth, snow coverage, and grain size).

Despite these advancements, limitations persist. The model relies on a simulated urban canyon structure for albedo analysis, suitable for studying snow effects but challenging for capturing the exact angular distribution for a specific district. The model simplifies the surface composition by assuming fixed snow proportions and materials, which may vary with factors like season, location, and study area. Accurate BRDF reproduction requires considering diverse surface characteristics. Additionally, the model assumes equal snow coverage on roofs and ground surfaces, which may not hold in real-world scenarios due to roof shape, slope, and snow removal. While initial validation with flat surfaces and MODIS data suggests the model’s effectiveness, a more comprehensive evaluation is hampered by limited multi-angle reflectance data. Future work will involve applying the model to a real urban structure with snow, using UVA technology for land surface parameter estimation and multi-angle observations to enhance model verification and applicability.

6. Conclusions

This study introduces a novel snow-covered urban canyon bidirectional reflectance model that simulates the impact of urban structures (urban coverage and building height), snow properties (snow coverage, snow grain size, and soot pollutant concentration), and solar zenith angles on the reflectance spectrum and angular distributions across seven MODIS bands (400–2400 nm). We validate the model using three reference cases and real scenarios from MODIS Terra and Aqua data in Harbin, demonstrating its ability to accurately reproduce spectral reflectance in snow-free, fresh snow, and snow melt conditions. A deeper analysis reveals that snow coverage is the primary determinant of albedo spectrum and bidirectional reflectance. Urban coverage, building height, and pollutant concentration have the most significant effect on visible and infrared band reflectance, while snow grain size influences shortwave infrared bands. As building height increases, bidirectional reflectance in the ALM shows less variation at relative azimuth angles between 120° and 180°. Snow coverage changes have minimal impact on the angular variation pattern of directional reflectance in the ALM, PP, and CP. However, at Band 6, bidirectional reflectance at low solar zenith angles (less than 30°) in the PP and CP is more susceptible to snow coverage variations. As snow grain size grows, the bidirectional reflectance at the backward scattering direction in Band 1 (0.5–0.6) is lower than the forward scattering direction (around 0.8) in the PP.

There are several limitations that require attention, including the use of an idealized city model, a fixed ratio of snow and man-made materials for reconstructing surface reflectance, and the assumption of uniform snow coverage on both roofs and the ground. In the future, we plan to apply our numerical model to a realistic 3D urban structure, incorporating more accurate surface reflectance and snow cover configurations to simulate the bidirectional reflectance distribution. This research is anticipated to contribute significantly to the understanding of how snow cover affects the anisotropic reflection characteristics of urban canyons. The outcomes may also catalyze further investigations in associated domains, such as quantifying the radiation interactions between snow and urban surfaces, assessing the sensitivity of urban reflectance to snow and building characteristics, differentiating between clouds and snow, and extracting aerosol optical depth along with snow parameters like coverage, grain size, and pollutant concentration from satellite observations.